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ABSTRACT  OF  THE  DISSERTATION 


Advances  in  Filtering  Techniques 
for  Stochastic  Systems  with  Uncertain  Parameters 

by 

Calvin  Chris  Schneider,  Jr. 

Doctor  of  Philosophy  in  Engineering 
University  of  California,  Los  Angeles,  1978 
Professor  Cornelius  T.  Leondes,  Chairman 


Optimum  estimation  of  the  states  of  linear  dynamic 
processes  using  noise  corrupted  measurements,  commonly  re- 
ferred to  as  Kalman  filtering,  requires  an  exact  knowledge 
of  the  equations  which  govern  the  complete  stochastic 
system.  In  practical  applications,  however,  these  equa- 
tions depend  on  parameters  which  can  not  be  precisely 
defined.  When  the  parameter  uncertainties  are  sufficiently 
large,  standard  filtering  techniques  can  produce  inaccurate 
and  inadequate  estimates.  In  this  dissertation,  two  alter- 
native techniques  of  state  estimation  for  systems  with  un- 
certain parameters  are  considered. 
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The  first  technique  assumes  a modification  to  the 
Kalman  filter  estimate  which  incorporates  the  parameter 


uncertainty  without  estimating  the  parameter.  The  second 
technique  defines  am  adaptive  estimate  of  the  uncertain 
parameters  which  is  then  utilized  to  improve  the  state 
estimates.  Both  techniques  are  applied  to  a numerical 
example.  Computational  results  are  generated  and  compared 
to  estimates  produced  by  a Kalman  filter  and  an  extended 
Kalman  filter.  The  effects  of  the  measurement  noise  level, 
a gap  in  the  measurement  data,  and  an  erroneous  data  spike 
on  the  filter  estimates  are  also  considered. 
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CHAPTER  1.0 


INTRODUCTION 

1.1  LINEAR  FILTERING 

Linear  filtering  is  concerned  with  finding  an  opti- 
mal estimate  of  some  quantity  (an  unknown  parameter,  a 
random  variable,  or  a random  signal)  when  an  adequate 
system  model  is  a linear  function  of  that  quantity.  The 
linear  equation  may  be  corrupted  by  additive  noise  and/or 
driven  by  known  inputs. 

One  of  the  first  studies  dealing  with  linear  esti- 
mation was  performed  in  the  early  1800's  by  Gauss  (1963) 
who  developed  the  method  of  least-squares  to  estimate  un- 
known parameters.  Legendre  is  also  considered  as  an 
original  inventor  of  the  same  theory  based  on  his  work  in 
1806  titled,  "Nouvelles  Methodes  pour  la  Determination  des 
Orbites  des  Cometes,"  (see  Sorenson  1970).  The  motivation 
in  the  development  of  estimation  theory  during  this  period 
was  based  on  the  desire  to  calculate  orbital  parameters  of 
the  motions  of  heavenly  bodies.  In  these  astronomical 
studies,  the  orbiting  bodies  were  observed  using  telescopes 
and  values  of  the  unknown  parameters  were  inferred  from 
the  measurement  data. 
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The  next  significant  development  in  filtering  theory 
was  achieved  in  the  1940' s by  Kolmogorov  (1941)  and  Wiener 
(1949)  for  the  class  of  problems  in  which  the  signal  and 
measurements  were  scalar,  random  processes.  Kolmogorov 
studied  only  discrete-time  problems  of  stationary  random 
processes  and  solved  them  by  using  a simple  representation 
of  such  processes  that  was  suggested  in  a doctoral  disserta- 
tion by  Wold  (1938) . This  representation,  which  is  obtained 
by  a recursive  orthonormalization  procedure,  is  known  as 
the  Wold  decomposition.  Wiener,  on  the  other  hand,  studied 
mainly  continuous-time  problems  and  took  an  almost  non- 
probabilistic  approach.  The  main  result  of  his  work  was  an 
integral  equation  called  the  Wiener-Hopf  equation.  The 
solution  of  this  equation  is  a weighting  function  which, 
when  convolved  with  the  corrupted  linear  measurement,  pro- 
duces an  unbiased  minimum  variance  estimate  of  the  random 
\ 

signal.  THe  Wiener-Hopf  equation,  however,  can  only  be 

* 

solved  explicitly  for  certain  stationary  random  processes. 
For  this  reason,  general  filter  design  using  the  Wiener- 
Hopf  integral  equation  has  only  limited  practical  applica- 
tions. Many  generalizations  of  this  work  were  pursued  in 
the  1940's  and  1950 's  by  Blackman,  Bode,  and  Shannon 
(1948)  , Bode  and  Shannon  (1950)  , Doob  (1953),  and  Lattin  and 
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Battin  (1958)  but  very  few  practical  results  were  achieved, 
notable  exceptions  being  the  work  of  Booton  (1952)  and 
Z6de}|  and  Ragazzini  (1950,  1952)  . 

In  order  to  circumvent  the  difficulty  of  solving 
the  integral  equation  for  the  general  filtering  problem, 
the  idea  of  generating  least-squares  estimates  in  a recur- 
sive manner  was  introduced  in  the  1950's.  This  interest 
was  stimulated  by  the  development  and  increased  use  of 
digital  computers.  With  the  digital  computer,  the  esti- 
mates were  generated  dynamically  by  a data  processing 
algorithm  which  represented  either  a differential  or  dif- 
ference equation.  The  first  works  of  recursive  least- 
squares  estimations  were  produced  by  Carlton  and  Follin 
(1956)  and  Swerling  (1959).  Swerling's  report  presented  a 
recursive  filtering  procedure  similar  to  that  described 
shortly  thereafter  by  Kalman  (1960)  and  is  generally  con- 
sidered to  be  the  starting  point  of  modern  estimation 
theory  commonly  referred  to  as  "Kalman  filtering." 

The  paper  by  Kalman  and  a subsequent  paper  by 
Kalman  and  Bucy  (1961)  generalized  the  results  of  linear 
estimation  theory.  The  main  contribution  of  these  papers 
was  the  transformation  of  the  Wiener-Hopf  integral  equation 
into  an  equivalent  nonlinear  differential  equation  whose 
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solution  yields  the  covariance  matrix  of  the  minimum  fil- 


tering error  and  provides  all  the  information  necessary  for 
the  design  of  the  optimal  filter.  Thus,  by  requiring  a 
numerical  rather  than  an  analytical  solution  to  the  Wiener- 
Hopf  equation,  a simple  recursive  filter  was  developed 
which  could  be  easily  synthesized  on  a digital  computer. 

The  optimality  of  the  Kalman  filter,  however,  is 
dependent  upon  complete  knowledge  of  the  parameters  which 
define  the  system.  In  practical  applications,  these  quan- 
tities may  not  be  precisely  defined.  When  the  parameter 
uncertainties  are  sufficiently  large,  the  standard  Kalman 
filter  can  produce  inaccurate  and  inadequate  estimates. 

1.2  SYSTEM  MODEL  AND  PROBLEM  STATEMENT 

Consider  a dynamic  system  represented  by  a linear 
difference  equation  and  driven  by  a combination  of  known 
inputs  and  white  noise.  The  system  may  be  a discrete-time 
process  or  a continuous-time  process  modeled  in  discrete 
time  by  means  of  a state  transition  matrix.  The  stochastic, 
vector  difference  equation  for  such  a system  is  given  by 

xi+l  = ♦i+l/i  (“}  Xi  + *i+l/i  (a)  ui 

+ (1.1) 
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where 


x ^ - n-dimensional  state  vector  at  time  instant  i 

- m-dimensional  input  vector  at  time  instant  i 

w^  - q-dimensional  noise  vector  at  time  instant  i 

- n x n transition  matrix  between  time 
instants  i and  i+1 

^x+l/i (a)  - n x m input  matrix  between  time 
instants  i and  i+1 

I i+l/i(<* ) - n x q noise  input  matrix  between  time 
instants  i and  i+1 

and  where  each  of  the  matrices  # ^i+l/i(a)»  and 

may  depend  on  p constant  parameters  designated  a. 
The  information  available  about  the  state  is  obtained  from 
measurements  related  to  the  state  according  to 

2 . = H . (a)  x . + v.  (1.2) 

1111 

where 

z - r-dimensional  measurement  vector  at 
time  instant  i 

v^  - r-dimensional  noise  vector  at 
time  instant  i 

- r x n measurement  matrix  at  time  instant  i 
and  where  the  measurement  matrix  H^(a)  may  also  depend  on 
the  parameters  a. 
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Although  the  same  parameters  may  not  exist  jointly 


X * 


in  the  measurement,  transition,  input,  and  noise  input 
matrices,  Hi(a),  (a) , and 

respectively,  this  representation  is  assumed  here  to 
generalize  the  results  to  follow.  The  case  involving 
different  parameters  in  each  of  the  matrices  is  then  a 
subset  of  this  general  representation.  Note  that  if  the 
actual  system  is  a continuous-time  system,  as  described  by 
the  equation 

x = Fx  + Bu  + Gw  (1.3) 


where 


F - continuous  state  transition  matrix 
B - continuous  input  matrix 
G - continuous  noise  input  matrix 

then  a single  uncertain  parameter  in  F will  generally 
affect  both  4>jy^_^(a)  and  ^i/i_i(a)  in  the  equivalent 
discrete  representation  (see  Appendix  A) . 

The  noise  vectors,  w^  and  v^,  are  assumed  to  be  un- 
correlated Gaussian  white  noise  sequences  with  statistics 


0 


(1.4) 


E[wJ  = 

E[vJ  = 0 

E Jw^Vj J = 0 for  all  i,j 
E[wiwj]  = QjAj 
E[vivi]  - Ri5ij 

where  *]•]  denotes  the  expectation  operator,  is  a 
positive  definite  q x q matrix,  is  a positive  definite 
r x r matrix,  and  8^^  is  the  Kronecker  delta  function. 

The  initial  state  xQ  is  defined  as  a Gaussian  ran- 
dom vector  defined  by  a mean  xQ  and  covariance  PQ  as 
follows: 


(1.5) 

(1.6) 

(1.7) 

(1.8) 


(1.9) 


e[(*o-Jo)(*<A)T]  - po  f1-10* 


where  P is  a positive  definite  n x n matrix, 
o 
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The  true  values  for  the  parameters  in  the  system 


r 

f 
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model  are  assumed  to  be  unknown,  the  expected  time- 
invariant  value  being  a0.  The  variation  Sa(i.e.,  a-aQ)  is 
assumed  to  be  a normally  distributed  vector  of  zero  mean 
and  variance  defined  by 

E £(8cr)  (5a)TJ  = cr  - diag  [o^,  o"2,  . . .crp]  (1.11) 

where  cr  is  a matrix  with  elements  along  the  diagonal 
designated  cr^  through  o-p. 

It  should  be  noted  that  parameter  uncertainties  in 
the  measurement  matrix  are  generally  small  or  nonexistent 
in  most  problems  compared  with  uncertainties  in  the  trans- 
ition and  input  matrices.  However,  they  are  assumed  here 
to  provide  a more  general  development. 

The  problem  is  to  obtain  consistent  and  improved 
estimates  of  the  states  x^  considering  the  uncertainty  of 
the  parameters  a,  the  objective  being  the  development  of 
advanced  techniques  of  state  estimation  which  avoid  the 
difficulties  of  previous  approaches. 

1.3  PREVIOUS  RESEARCH 

Originally,  the  problem  of  state  estimation  for 
stochastic  systems  with  uncertain  parameters  was  posed  as  a 
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non linear  estimation  problem  by  augmenting  the  state  vector 
with  the  unknown  parameters  in  the  state  equations.  The 
extended  Kalman-Bucy  filter  was  then  applied  to  the  re- 
sulting nonlinear  filtering  problem  by  Kopp  and  Orford 
(1963) , Farison,  Graham,  and  Shelton  (1967) , and  Jazwinski 
(1970) . Although  this  technique  is  feasible,  the  computa- 
tional burden  is  increased,  and  the  filter  is  subject  to 
divergence.  The  divergence  is  due  primarily  to  the  linear- 
ization about  the  current  estimate  of  the  extended  Kalman 
filter.  If  the  a priori  estimates  are  poor  or  if  the 
filtered  estimates  lie  outside  the  linear  region,  the 
filter  will  often  diverge. 

Another  approach  by  Neal  (1967)  has  been  to  arti- 
ficially increase  selected  elements  of  the  plant  noise 
covariance  matrix  to  reflect  the  increased  uncertainty  in 
the  state  dynamics.  However,  this  trick  is  a process  of 
judiciously  tuning  the  Kalman  filter  and  becomes  more  of  an 
art  than  a science.  Recently,  a trajectory  sensitivity 
approach  was  taken  by  Chang  and  Belanger  (1976)  to  the 
design  of  the  Kalman  filter  gain  for  a continuous-time 
system.  A two-point  boundary  value  problem  was  formulated 
where  the  performance  index  was  a quadratic  function  of  the 
final  error.  Through  expedient  approximations,  a solution 
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was  obtained  and  the  final  estimate  was  more  accurate  than 
that  of  the  nominal  Kalman  filter;  however,  intermediate 
estimates  did  not  show  as  much  improvement. 

In  adaptive  filtering,  the  problem  has  been  formu- 
lated using  decision  theory  by  Magill  (1965)  and  Lainiotis 
(1971)  and  recently  using  maximum  likelihood  by  Maybeck 
(1968),  Mehra  (1970),  and  Gupta  and  Mehra  (1974).  The 
decision  theory  approach  assumes  a bank  of  parallel  Kalman 
filters  with  each  filter  matched  to  a possible  parameter 
vector  value.  The  state  estimates  generated  by  the  Kalman 
filters  are  then  combined  using  a weighted  Siam  in  which  the 
Bayesian  a posteriori  hypothesis  probabilities  are  weighting 
factors.  If  one  of  the  selected  parameter  vectors  coin- 
cides  with  the  true  parameter,  this  method  gives  the  mini- 
mum variance  estimates  of  both  the  state  and  parameters. 

The  computational  burden  of  this  method,  however,  can  be 
quite  significant.  The  maximum  likelihood  estimation  is 
based  on  the  innovation  sequence  of  the  system.  However, 
since  the  parameters  and  the  parameter  probability  distri- 
bution are  both  assumed  to  be  unknown,  the  computation 
aspects  of  the  solution  can  be  extremely  difficult. 

1.4  PROPOSED  APPROACH 

The  research  presented  herein  consists  of  developing 
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advanced  techniques  of  state  estimation  for  stochastic 
systems  with  uncertain  parameters . Both  nonadaptive  and 
adaptive  techniques  are  considered.  The  nonadaptive 
approach  assumes  a Kalman-type  filter  defined  as 

*i  = *i  + Gi  | zi_Hi(ao)^i  | + A»i  (1.12) 


where 


xi  = <t>. 


i/i-l(ao>*i-l  + Vi-l'V"!-!  (1-13) 


and 


x^  - predicted  value  for  x^_ 

Gi  - modified  Kalman  gain  matrix 
- estimator  bias 

The  transition  matrices  of  the  state  equation  are  expanded 
to  first  order  in  the  a parameters.  The  state  equation 
becomes 

xi  * xi(a0)  + 


P dx . 


:=1 


das 


8ai 


(1.14) 


11 


where 


xi(“o)  * *i/i-l(oo)xi-l  + '"i/i-l(“o>ui.l  + ri/i-l<“o>wi-l 

(1.15) 


da. 


d*i/i-l  . dVi-l  d/i/i-l 

~"d^  Xi-1  + ~ Ui-1  + Wi-1 


lak 


lak 


(1.16) 


The  bias  of  the  estimator  is  then  determined  by  requiring 


that  the  expected  value  of  the  error  e^  (i.e.,  x^-x^)  be 
zero  or 


f*i  = f(E[eJ  = 0) 


(1.17) 


The  filter  gain  is  defined  by  requiring  that  the  expected 
value  of  a quadratic  function  of  the  estimation  error  is 
minimized.  Equivalently,  this  is  given  by 


In  the  adaptive  technique,  a Kalman  filter  is 
coupled  with  a Bayesian  estimator  of  the  uncertain 
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parameters.  The  Bayesian  estimator  is  a maximum  a 
posteriori  estimator  obtained  by  maximizing  the  conditional 
density  function  of  the  parameters  given  by 


a 


= f 


op  (<*|z) 


8a 


= 0 


(1.19) 


where  p(a|z)  is  the  conditional  density  function  of  the 
parameters  given  the  measurements  Z.  The  state  estimate  is 
then  determined  using  the  estimate  of  the  parameters.  The 
state  estimate  is  given  by  the  standard  Kalman  filter 
equation 


(1.20) 


where  Ki  is  the  Kalman  gain,  and  where  all  the  terms  that 
contain  the  parameters  a are  determined  using  the 
estimate  a. 

1.5  SUMMARY  OF  DISSERTATION 

Chapter  2 develops  the  Kalman  filter  algorithm  for 
linear,  discrete  processes.  The  problem  of  uncertain  para- 
meters in  the  system  equation  is  then  posed  as  a nonlinear 
problem  by  adjoining  the  uncertain  parameters  to  the  state 
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vector,  and  the  extended  Kalman  filter  is  derived  as  a 


solution  to  this  nonlinear  problem. 

In  Chapters  3 and  4,  research  on  alternative 
methods  of  state  estimation  for  stochastic  systems  is  pre- 
sented. Chapter  3 provides  the  derivation  of  a modified 
Kalman  filter,  and  complete  estimation  algorithms  are 
developed.  A discussion  of  the  stability  of  the  modified 
filter  is  also  provided.  Chapter  4 develops  an  adaptive 
estimator  of  the  parameters,  which  is  then  utilized  to 
improve  the  state  estimates.  In  addition,  the  asymptotic 
nature  of  the  parameter  estimator  is  assessed. 

The  techniques  presented  in  this  dissertation  are 
applied  to  a practical  problem  in  Chapter  5.  Computational 
results  are  generated  for  the  filters  proposed  in  Chapters 
3 and  4.  A comparison  of  these  results  with  those  produced 
by  the  Kalman  filter  and  the  extended  Kalman  filter  in 
Chapter  2 is  also  presented.  Finally,  Chapter  6 presents 
the  conclusions  of  the  research  and  suggestions  for  future 
investigations . 
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CHAPTER  2 

KALMAN  FILTER  AND  EXTENDED  KALMAN  FILTER 

This  chapter  develops  the  Kalman  and  extended 
Kalman  filters  for  solving  the  problem  of  state  estimation 
from  samples  of  noise-corrupted  data.  The  development  is 
based  on  the  assumption  that  all  parameters  of  the  system 
structure  are  known  precisely.  Section  2.1  derives  the 
Kalman  filter  for  systems  represented  by  linear  difference 
equations.  The  maximum  likelihood  concept  is  utilized  to 
define  the  likelihood  function  from  the  available  statis- 
tical information  and  this  function  generates  the  likelihood 
equation.  The  solution  of  the  likelihood  equation  provides 
the  desired  Kalman  filter  expression.  In  Section  2.2,  the 
problem  of  uncertain  parameters  in  the  system  equations  is 
considered.  The  linear  problem  is  transformed  to  a non- 
linear problem  by  augmenting  the  uncertain  parameters  to 
the  state  vector.  The  extended  Kalman  filter  is  then 
derived  to  solve  this  nonlinear  problem. 

2.1  OPTIMUM  STATE  ESTIMATION 

This  section  develops  the  equations  of  the  Kalman 
filter,  the  statistics  of  the  filter  estimation  error,  the 
statistics  of  the  measurement  information,  and  the 


I 
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statistics  of  the  measurement  residual.  The  system  model 


f 


! 

[ 


used  for  this  development  was  presented  in  the  previous 
chapter  with  the  assumption  that  the  true  value  of  the 
parameters  in  the  model  are  known. 

2.1.1  Derivation  of  the  State  Estimate 

Let  Z^  denote  the  vector  of  realized  values  of  the 

measurements  z,  ...  z.;  thus,  Z.  is  a vector  of  dimension 
1 i i 

(i*m) . In  order*  to  obtain  a maximum  likelihood  estimate  of 
the  system  state,  a likelihood  function  must  be  defined 
relating  the  measurements  (whose  values  are  known)  and  the 
state  variables  (the  unknowns  to  be  estimated) . An  appro- 
priate likelihood  function  is  the  conditional  probability 
density  of  the  state  given  the  measurements  p(x^|z^). 

(Note  that  all  system  parameters  are  implicitly  assumed  to 
be  known.) 

To  specify  the  maximum  likelihood  estimate,  the 
value  x^  is  sought  that  will  yield  the  greatest  possible 
value  of  the  likelihood  function.  In  this  particular  case, 
x^  is  the  value  that  will  yield  the  peak  of  the  conditional 
density  function  p(x^|z^),  i.e.,  it  will  maximize  the  pro- 
bability of  the  variables  to  be  estimated,  conditioned 
upon  the  events  actually  observed.  Rather  than  directly 
maximize  p(xjz^),  it  is  more  convenient  to  find  the 


16 


1 


minimum  of  the  negative  natural  logorithm  of  p(xi,Z^). 
This  is  a valid  step  because  for  any  function  p,  p and 
-ln(p)  attain  their  maxima  and  minima,  respectively,  at 
the  same  point.  Thus,  the  log- likelihood  function  is 
defined  as 


L(xi,Zi)  = -In  p(xi|zi)  (2.1) 


in  order  to  evaluate  the  desired  state  estimate. 

Using  Bayes'  rule,  the  conditional  density  is  first 
written  into  a more  convenient  form: 


p(x.,Z.) 

i i 

P(Zi> 


p(2.|Zi_i)  ptZ^) 


p(2i|Zi_1) 

By  means  of  this  equation,  the  conditional  density  for  x^, 
given  all  measurements  through  z^,  is  directly  related  to 
the  density  conditioned  on  all  measurements  up  to,  but  not 


(xj  Zi-1)  (2*2) 
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including,  z^.  It  is  used  to  relate  the  estimate  of  x^ 
just  after  the  measurement  to  the  estimate  just  before 
z^.  The  component  terms  of  the  above  equation,  however, 
must  first  be  evaluated. 

Let  the  mean  of  the  density  p(x^|z^  be  denoted 
as  x^  and  the  corresponding  conditional  covariance  as 
pi/i_l*  Then, 


‘i  - 


(2.3) 


•i-i = E[(*i-Xi>  (xi-*i>Tlzi-i] 


(2.4) 


These  expectations  are  ensemble  averages  over  all 
possible  initial  conditions  and  both  driving  and  measurement 
noise  sequences,  all  conditioned  on  the  values  of  the  mea- 
surements. Because  the  system  is  linear,  driven  by  a white 
Gaussian  noise  sequence  and/or  a deterministic  input,  and 
the  probability  density  function  of  the  initial  state  is 
Gaussian,  the  conditional  probability  density  function  of 
x^  before  the  measurement  at  instant  i is  normal.  Thus, 
the  density  function  is  given  by  the  equation 
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exp 


(2^)  I pi/i_1l 


-*S  Cxi-x  J Tpi/i-i  [xi-xi] 


(2.5) 


where  J I is  the  detern,inant  of 

Consider  now  the  term  pfz^fx^.Z^  The  observation 

at  time  i is 

z . = H.x.  + v.  (1.2) 

1 1 i x 

(Note  that  the  term  (a)  is  shortened  to  since  the 

parameters  are  assumed  to  be  known.)  The  sequence  has 

been  defined  to  be  normally  distributed  and  is  independent 
of  x^.  Therefore,  conditioned  on  a particular  value  of  x^, 
is  a normal  random  variable  with  conditional  mean 

E[*Jvzi-i]  * 4iK] 

I - HiXi  (2.6) 

and  conditional  cbvariance 
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E[l>i-ViH*i-Hi*i]T|*i.Zi_1]  * e[vavT]  = Ra  (2.7) 


so  that  the  conditional  probability  density  of  z^,  given  x^ 
and  Zi_^,  is 


P<zilxi'zi-1>  * — 7R  exp 

(2tt)  iRil 


= -Jsfz.  -H.x.]TR71f'z.  -H.x.]  | (2.8) 

LiiiJ  1 L i 1 iJ  I 


Finally,  consider  p(z.JZ^  ^)  . The  output  of  the 
dynamic  equation  is  a Gauss -Markov  sequence;  in  other  words, 
x^  is  a Gaussian  random  variable.  Furthermore,  linear 
operations  on  Gaussian  random  variables  produce  Gaussian 
random  variables,  so  H^x^  is  normally  distributed.  The  sum 
of  two  independent  normal  random  variables  (i.e.,  H^x^  and 
v^)  yields  a normal  variable.  Therefore,  p(zi|z^_j_)  is  a 
normal  density  characterized  by  conditional  mean  and 
covariance : 


E [zii zi_i]  - HtEh|zi-i] + E[vi|zi.i]  ■ Vi  <2-9> 

E [Wi*  (zi-Hi*i>Tlzi-i]  ■ Hipi/i-iHI + Ri  <2-10> 


and  given  by  the  equation 


P<Hlzi-i> J75 1 


<2*>  lHipi/i-iHi  + Ril 


* |-)S<zi-Hi*i>T<HiP1/i_1H?  + R^-1  ' 


(2.11) 


Parenthetically,  note  that  this  term  is  not  a function  of  x 
and  therefore  will  have  no  effect  on  the  derivation. 

By  incorporating  the  above  equations  into  the  ex- 
pression for  p(xjjz^)  and  taking  the  natural  logarithm,  the 
negative  log- likelihood  function  is  obtained  as 


MXi.Zi)  = }5/n|(27r)n|HiPi/i_1Hi  + RiT^R*)  I Pi/i.J  | 


+ ?5  J(x.-x.  )TP.  . (x.-x.) 
•'l  l 1 l/l-l  1 l' 


+ (zi-Hixi)TRT1(zi-Hixi) 


ij|(Zi-HiXi)T(HiPi/i_iH^  + Ri)  ”1(zi-Hixi) 
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where  the  first  and  last  terms  do  not  involve  x^. 

To  obtain  the  maximum  likelihood  estimate  of  x(i), 
the  equation 


dL(x.  ,Z.) 

i i 


is  solved,  where  x^  denotes  the  maximum  likelihood  estimate 
of  x^.  Substituting  L(x^,Z^)  into  this  equation  yields 


xf 


A 

*-X  . 


= 0 


(2.13) 


pI/i-i(xi-xt) 


htrT1(z.-h.x. ) 
11  1 11 


(2.14) 


The  solution  is  then 


+ 'l  <PI/i-lXi  + <2'15> 


Instead  of  solving  the  likelihood  Equation  (2.13) 
directly  by  using  Equation  (2.12),  algebraic  minipulations 
(completing  squares  and  applying  the  Householder  inversion 
lemma)  can  first  be  performed  upon  the  negative  log- 
likelihood  function  to  obtain  the  form 
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L [x  (i)  ,Z  (i)]  = 35/n  | (2tt) nJ  P±|  | 


+ H(x.-x.)Tp:1(x.-x.) 
i 1.  i.  i i 


where 


ft.  = 


(ptV  , + htr"1h.)"1  • (PT y.  ,X.  + htr-1z.) 

i/i-l  111  1/1-I  1 111 


pT1  = p7/.  + htr71h. 

1 1/1-1  111 


(2.16) 


(2.17) 


(2.18) 


This  demonstrates  that  pfxJZjJ  is  a normal  density  with 


mean  x^  and  covariance  P^: 


p(xi|zi)  = 


n/2 , , h 

(2tt)  |pil 


exp  | • 


= -J5(xi-xi)TP~1(xi-£i)  (2.19) 


ft.  = E fx . | Z . ] 
1 L 1 1 


Pi 


= Ei[(xi-fti)  (xi-fti)T] 


(2.20) 


(2.21) 
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As  expected,  due  to  the  symmetry  of  this  density,  the 
maximum  likelihood  estimate  and  the  mean  of  the  density 
coincide.  The  matrix  inversion  lemma  can  be  used  to  put 
this  estimate  into  a better  computational  form,  requiring 
inversion  of  an  m x m matrix  instead  of  an  n x n (m  is 
typically  smaller  than  n and,  since  measurements  can  be 
incorporated  one  at  a time,  the  inversion  can  become  a 
simple  division).  F.or  positive  definite  and  , the 

lemma  yields 


(PTy.  . + htr-1h  . ) -1  = p.  ..  , 
1/1-1  iii  1/1-1 


(2.22) 


so  that  the  state  estimate  can  be  written  as 


xi  + Ki(zi~Hixi> 


(1-20) 


where 


K. 

i 


= Pi/i-lHi + Xi) 


-1 


(2.23) 
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which  can  be  recognized  as  the  form  of  the  Kalman  filter 


with  K^  as  the  Kalman  gain.  The  update  of  the  covariance 
matrix  from  just  before  the  i^  measurement  to  just  after 
is  given  by  Equation  (2.18)  rewritten  as 


P. 

i 


(P7/.  + hTr-^.)”1 

i/l-l  ii  i 


(2.24) 


which,  although  it  adds  the  measurement  information  in  a 
simple  way,  requires  an  n x n matrix  inversion.  Equation 
(2.22)  yields 


pi  3 Vi-i  - pi/i-iHi  + Ri'~lHipi/i-i 


(2.25) 


This  expression  can  also  be  written 


■ I'-WVu 


(2.26) 


Although  these  forms  involve  m x m inversions,  they  do  have 
some  undesirable  characteristics.  Equation  (2.25)  is  some- 
times a small  difference  of  large  numbers  (especially  if  the 
measurements  are  very  accurate)  and  thus  is  subject  to 
numerical  precision  problems.  Equation  (2.26)  does  not 
assure  positive  definiteness  or  symmetry,  and  thus  can  also 


-J 


lead  to  numerical  difficulties.  Now,  if  the  estimate  in 
Equation  (1.20)  is  rewritten  as 


A 

x. 


= (i-K.H.)x. 

ill 


K.z. 

l l 


(2.27) 


then  it  can  readily  be  shown  that  an  equivalent  expression 
for  would  be 


P. 

i 


(I-K.H.)P./.  . (I-K.H.  )T  + K.R.kT 

l l i/l-l  ll  ill 


(2.28) 


This  is  the  sum  of  two  symmetric,  positive  definite  matrices 
so  numerical  computations  based  upon  this  form  will  be 
better  conditioned,  better  assuring  the  symmetry  and 
positive  definiteness  of  P^. 

The  probability  densities  involved  are  normal,  and 
the  propagation  of  the  entire  density  function  can  be 
specified  by  the  time  history  of  its  mean  and  covariance. 
Thus,  estimates  propagated  between  measurements  are 


Xi+1  - ECxi+ilzJ 


(2.29) 


or 
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(2.30) 


and 


xi-i  ■ WA  + ,"i+i/iui 


Pi+l/i  E [*xi+l_xi+l*  *xi+l~xi+l>T|zi] 

' E[k+l/i(xiA>  + ri+l/iwi|- 

•|Wi(vA  + ^i/i“i|T|zi]  <2 

or 

Pi+l/i  = ^i+l/ipi^i+l/i  + /i+l/iQiri+l/i  (2 

where  ^i+l/i(a)'  and  ^i+l/i^a^  are  written 

^i+l/i»  ^i+i/i*  and  ^2+1/i'  respectively,  since  the 
parameters  a are  assumed  to  be  known. 

To  summarize,  at  a measurement,  the  estimate  is 
updated  using; 


xi  - + K.z. 


(2 


P. 

i 


= (I-K.H  . ) P . ..  (I-K.H.)T  + K.R.kT 

i i i/i-l  ii  ill 
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(2 


.31) 

.32) 


27) 

28) 


A — — — 


Main  i 


where 


K. 

i 


T T -1 

P.  .H1  (H.P.  HT  + R . ) 
l/l-l  l l l/l-l  l l 


(2.23) 


and,  to  propagate  the  estimate  up  to  the  time  of  the  next 
measurement,  the  relations  are 


X.  -i  = <t>-  . + 'll.  , /.U. 

l+l  i+l/i  i i+1/i  i 


(2.30) 


Pi+l/i  ^i+l/iPi^i+l/i  + ri+l/iQiri+l/i  (2'32) 


Note  that  can  also  be  expressed  as 


K.  = P.HTR-1 

l ill 


(2.33) 


(This  can  be  expanded  using  the  expression  for  P^  to  show 
it  is  equal  to  the  previous  evaluation  of  K^.)  Although 
it  is  a simpler  expression  and  has  the  same  appearance  as 
the  continuous -time  Kalman  gain,  it  is  not  as  convenient  to 
use  in  recursions  (see  Ho  1962) . 

2.1.2  Statistics  of  the  Estimation  Error 

This  section  defines  the  statistics  of  the  error 
committed  by  the  estimate  immediately  after  the  measurement 


X * 


at  time  i.  The  error  is  given  by  the  equation 


I 


I 

A 
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It  should  be  recognized  also  that,  since  the  expected  value 
of  e^  is  zero,  the  state  estimate  x^  is  unbiased. 

2.1.3  Statistics  of  the  Measurement  Information 

Consider  the  statistics  of  the  new  state  informa- 
tion introduced  at  each  measurement.  From  Equation  (1.20), 
this  information  is  portrayed  by 


x.  -x.  =K.(z.-H.x.) 
1 1 1111 


(2.38) 


Upon  substitution  of  z^,  this  equation  becomes 


Wi-Vi + vi> 


+ + Vi 


(2.39) 


where  the  term  in  brackets  is  the  sum  of  three  independent 
terms,  each  independent  of  ..  Consequently,  the  value 
(A^-x^)  is  independent  of  Thus,  the  statistics  of 

state  information  added  by  a measurement  are  expressed  as: 
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E[<*i-*i>  2i-l]  * E(vxi] 


(2.40) 


E [(Vxi>  <xi-xi)T|Zi-l]  * E [(xi-xi>  (ii-xi)T] 


Ki<HiPi/i-lHi  + Ri>KI 


= P...  HT  (H , P . HT  + R.)-1H.P.  ..  . 
1/1-1  1 1 l/l-l  111  l/l-l 


= K.H.P.  ..  . 
l l l/l-l 


(2.41) 


Recognize  that  the  state  information  covariance  is  the 
quantity  that  reduces  the  state  estimate  error  covariance 
at  each  measurement  (see  Equation  2.26) 

2.1.4  Statistics  of  the  Measurement  Residual 

The  measurement  residual  at  time  is  is  the  differ- 
ence between  the  measurement  and  the  best  prediction  of 
zi  based  on  Z^.^,  Hi*i*  This  quantity,  given  by  the 
equation 


"i  - zi  - Hixi 


(2.42) 


1 


may  be  regarded  as  the  new  information  from  the  current 

observation  z.,  given  all  the  past  observations  Z.  , and 
i i-l 

the  past  information  deduced  from  them.  As  such,  the  term 
is  commonly  referred  to  as  the  innovation  and  the  series 
of  residuals  is  defined  as  the  innovation  sequence.  This 
designation  was  first  used  by  Wiener  to  describe  such 
processes . 

The  innovation  represents  a linear  operation  on 
Gaussian  random  variables  and,  therefore,  is  itself  a 
Gaussian  random  variable.  The  mean  ~v.  and  variance  of  the 

i 

innovation  V\  are  given  by 

V i = E(zi-Hixi|zi_1) 

= E(HixijZi_1)  + E(vi)Zi_1)  - E(Hixi|zi_1) 

= 

= 0 (2.43) 
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follows : 


E[€.€T  .] 

L i 1-1J 


can  be  evaluated  as 


For  i not  equal  to  i-j,  the  term  E [v^vT  is  zero. 
Therefore,  the  correlation  function  becomes 


* 0 


(2.48) 


(2.49, 


This  equation  demonstrates  that  the  innovation  process  is 
a white-noise  process  (i.e.,  the  innovations  are 


independent  from  time  point  to  time  point) . This  outcome 
is  the  result  to  be  expected  for  an  optimum  filter. 


It  is  interesting  to  note  that,  while  the  innova- 
tions are  statistically  independent  as  shown  above,  they 
are  related  to  the  error  of  the  predicted  state  error  e ^ 
which  can  be  defined  by  a recursive  relationship.  The 
equation  relating  the  innovation  and  predicted  state 
error  is 

^i  = zi  — Hixi 


= H. (x.-x. ) + v. 

111  i 

= H.  €.  + v.  (2.50) 

11  i 

The  predicted  state  error  can  be  written 

e.  = x.-  x. 

i ii 

* Vi-lxi-l  + ri/i-lwi-l  - ♦i/i-A-l  <2'51) 

Substitution  for  x^  ^ gives  the  desired  recursive  relation- 
ship 
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€i  " ^i/i-l  [1-Ki-lHi-l]  (xi-l”xi-l) 

+ ri/i-lwi-l  ~ Vi-lKi-lvi-l 

= *i/i-l[I_Ki-lHi-l]  €i-l 


+ -w.  , - <t>.  / . . K . .v.  , 

1 i/i-l  i-l  i/i-l  i-l  i-l 


(2.52) 


2.2  EXTENSION  OF  THE  OPTIMUM  STATE  ESTIMATE 

This  section  transforms  the  linear  difference 
equations  of  a stochastic  system  with  uncertain  parameters 
to  nonlinear  difference  equations.  The  extended  Kalman 
filter  is  then  derived  to  solve  this  nonlinear  estimation 
problem. 

2.2.1  Nonlinear  System  Equations 

In  the  preceding  section,  the  Kalman  filter  was 
derived  for  a stochastic  system  represented  by  linear 
difference  equations.  The  derivation  was  based  on  the 
assumption  that  the  parameters  in  the  system  equations  were 
known.  If  these  parameters  are  unknown,  they  may  be 
treated  as  additional  states  to  be  estimated;  however,  the 
system  equations  become  nonlinear. 
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Consider  augmenting  the  state  vector  with  the 
unknown  parameters  to  form  a new  state  vector  x|  given  by 


(2.53) 


The  state  and  measurements  may  then  be  rewritten  such  that 


xi  = f(xi-l'ui-l>  + 9(xi-l)wl-] 


Zi  = h(xj.)  + vi 


(2.54) 


(2.55) 


where  i + ^i/i  l^“^ui  1 ^ s inherent  in 

f(x|  ^.u^  and  H i(a)x^  is  represented  by  h(x|).  'These 
equations  are,  in  general,  nonlinear,  and  thus  provide  the 
basis  for  the  derivation  of  the  extended  Kalman  filter. 
2.2.2  State  Estimation  for  Nonlinear  System  Equations 


Assume  that  the  measurements  ^ are  available  so 
that  the  optimum  state  estimate  at  time  i-1,  x^  and  the 
predicted  state  estimate  at  time  i,  x|,  are  provided.  The 
state  and  measurement  equations  can  then  be  expanded  in  a 
Taylor  series  about  these  estimates.  Introducing  the 
expansions  truncated  to  zero  or  first  order. 
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f(x!  ,u.  .)  = f(x:  ,u.  .)  + ~~ — 

1-1  1-1  1-1  1-1  aA,^ 


h-i-Si-i]  <2-56> 


g(xU)  - g(ftU) 


(2.57) 


h (x ! ) =h(x:)  + SiL.  (x!  -x! ) 

i i ax|  1 1 


(2.58) 


yields  new  system  equations: 


4>\  ..  ,x!  ..  + u!  . + T'  /■  iw-  i 
i/i-l  i-l  i-l  1 i/i-l  l-l 


(2.59) 


where 


z ! = H ! x!  + v.  , 

i ii  i-l 


(2.60) 


<b> 

-i/i-l 


5f(x:  ,u.  ) 

i-i  i-i 

dx‘ 

oxi-l 


(2.61) 


df(x.  wU.  .) 

ui-i  = f (fti-i'ui-i) — — ^i-i 

dxi-i 


(2.62) 


ri/i-i  * 9(XU> 


(2.63) 


h(x' ) 


(2.64) 


+ 


ah(x')  _ 

= — 

ax|  1 


ah(x;) 

i 

dx\ 

i 


(2.65) 


By  comparison,  these  system  equations  are  identical  in 
form  to  the  basic  equations  used  to  derive  the  Kalman 
filter  (i.e.,  the  equations  are  linear  with  known  matrices 
and  known  inputs) . Therefore,  the  estimate  of  x|  can  be 
developed  analogously  using  the  Kalman  filter  algorithms. 
The  estimate  of  x|  after  i measurements  is  given  by  the 
equation 


a, 

x! 

i 


x:  + k: (z: 

i i i 


H ! x ! ) 

i i 


X!  + K! 


Z • - 


dh(x!) 

h(x!)  + 

d*i 


1 X! 
X1 


dh(x!)  _ 


dxl 


*i 


- *i + Ki[*i  - h(xi>] 


(2.66) 


where  the  predicted  value  of  x|,  given  Z ^ ± measurements, 
is 


= f(x!_1,Ui_1) 


(2.67) 
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Equations  (2.66)  and  (2.70)  provide  the  complete  set  of 
relationships  for  the  extended  Kalman  filter.  Estimation 
of  both  the  state  x^  and  the  parameters  a is  accomplished 


by  this  filter. 
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CHAPTER  3 


MODIFIED  KALMAN  FILTER 
FOR  SYSTEMS  WITH  UNCERTAIN  PARAMETERS 

This  chapter  proposes  a modification  to  the  Kalman 
filter  for  state  estimation  of  discrete-time  systems  when 
the  parameters  in  the  system  equations  are  not  known 
exactly.  The  modification  includes  the  addition  of  a bias 
term  to  the  filter  and  a change  in  the  filter  gain  to 
reflect  the  uncertainty  in  the  system  parameters.  The 
error  covariance  matrices  of  the  filter  are  also  amended 
as  a result  of  the  parameter  uncertainties.  The  modifica- 
tion to  the  Kalman  filter  is  based  on  the  assumption  that 
the  parameter  uncertainties  are  relatively  small  so  that 
the  terms  containing  the  uncertain  parameters  can  be  ex- 
panded to  first  order  in  these  parameters.  Section  3.1 
delineates  this  expansion  and  provides  the  equational  rela- 
tionships that  are  utilized  in  Section  3'. 2 to  define  the 
modified  Kalman  filter. 

3.1  SYSTEM  MODEL  MODIFICATION 

Parenthetically,  it  should  be  recognized  that  a 
fundamental  assumption  in  this  development  is  that  a linear 
model  has  been  defined  which  adequately  represents  the 
behavior  of  a physical  system,  but  that  certain  parameters 


- 

J 
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in  this  linear  description  are  not  know  exactly.  Thus, 
equations  of  the  model  depict  those  aspects  of  a physical 
process  that  are  most  significant  to  a user's  purposes. 

The  filter  consequently  does  not  try  to  estimate  the 
"true"  values  of  the  states  of  the  system,  but  attempts 
to  determine  the  values  which,  when  substituted  into  the 
model,  yield  a model  output  behavior  that  best  duplicates 
the  actual  system  performance  (which  can  be  measured  only 
with  limited  certainty) . Moreover,  this  optimal  replica- 
tion is  achieved  for  a single  set  of  input  and  output 
sequences,  with  no  guarantee  that  the  same  model  provides 
an  optimal  system  representation  for  other  sequences  as 
well. 

The  ability  to  perform  the  estimation  will  thus 
be  determined  by  the  conditions  by  which  the  mathematical 
model  serves  as  an  adequate  system  representation.  If  the 
estimation  is  impossible  with  the  originally  formulated 
problem,  one  may  need  to  incorporate  different  measurements, 
additional  measurements,  or  alter  the  system  model,  in 
order  to  satisfy  these  conditions. 


X / 


However,  given  that  system  model  stipulated  in 


Chapter  1 is  adequate  and  that  the  variances  of  the 
parameters  are  sufficiently  small,  the  matrices 

containing  the  uncertain  parameters  can  be  expanded  to 
first  order  about  the  nominal  value  of  the  parameters 
designated  generically  aQ.  The  transition,  input,  and 
noise  input  matrices  become 


Vi-iw  - Vi-i'V  *LA/iA  (3-l> 

k=l  K 

*W ■ Vi-l'V  +E  «v P-2) 

k=l  K 

+E„/i/i-i5“k  (3-3> 

k«l  k 


where  the  lower  left  subscript  a indicates  differentiation 

1c 

with  respect  to  the  parameter.  Substituting  into  the 
state  equation  gives 


(3.4) 


+ [ri/i-l  (“o)  +£  akr L/i-l6ak]wi-l 


k=l 


Equation  (3.4)  may  be  separated  into  two  components 


xi  “ xi  <%>  +S 


(3.5) 


k=l 


where 


x . (a  ) = <t>.  . . (a  ) x . + ifi.  ..  . (a  ) u . + H . (a  ) w . 

l o l/l-l  o i-l  i/i-l  o i-l  i/i-l  o i-l 


(3.6) 


x.  = <6.  , x.  , + (/».,.  .u  i + /"•  /.  .w.  . (3.7) 

a.  l i/i-l  i-l  a./  i/i-l  i-l  a.  i/i-l  i-l 


and  where,  again,  the  lower  left  subscript  indicates  dif- 


.th 


ferentiation  with  respect  to  the  k parameter  The 


first  term  in  Equation  (3.5)  defines  the  transition  of  the 
states  and  the  effect  of  the  inputs  from  time  i-l  to  time  i 


based  on  nominal  values  for  the  parameters  a^.  The  second 


term  provides  the  change  in  the  states  at  time  i due  to 
the  variation  8 a from  the  nominal  value  a . 


V 


r 


The  measurement  matrix  can  also  be  expanded  to 
first  order  and  written 


Hi  (a)  = H.(«o)  +S<,vHiSak 

k k 


(3.8) 


where  the  lower  left  subscript  indicates  differentiation 
with  respect  to  the  kth  parameter.  The  measurement 
equation  is  thus  given  by 


:i  -[w  +L%His“kh + vi 


(3.9) 


and  utilizing  Equation  (3.4),  becomes 


!i  * Hi(“o)[',i/i-l(ao,xi-l  + *i/i-l(°o)Ui-l  + ri/l-l(ao)Wi-l] 

% V-otWi-l  + a + 

J\  K J* 

£ELHiU/i-l*i-l  + + a/i/i-lwi-l]  S“/S“k 

^ kit  “k  k k I 


+ vi 


(3.10) 
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expressions  for  the  system  model,  which  enable  modification 
of  the  Kalman  filter  and  permit  derivation  of  a state 
estimator  for  stochastic  systems  with  uncertain  parameters. 
3.2  DERIVATION  OF  THE  MODIFIED  KALMAN  FILTER 

The  Kalman  filter  is  an  optimum  estimator  for 
linear  stochastic  systems  only  when  the  parameters  of  the 
system  are  known.  When  the  system  parameters  deviate  from 
assigned  values,  the  decreased  accuracy  in  the  estimate 
computed  by  the  Kalman  filter  may  be  unacceptable.  If 
modifications  to  the  Kalman  filter  are  considered,  an 
improved  estimate  can  be  made  which  may  be  sufficiently 
accurate. 

How  shall  the  Kalman  filter  be  modified?  The  most 
likely  values  for  the  parameters  are  the  expected  values 
and,  thus,  terms  containing  these  parameters  should  be 
determined  with  these  nominal  values.  However,  since  the 
true  parameter  values  may  differ  from  the  nominal  values, 
the  modified  filter  may  be  biased.  Therefore,  a term  should 
be  added  to  the  filter  equation  to  compensate  for  this  bias. 
In  addition,  the  filter  gain  should  be  derived  to  reflect 
the  uncertainty  in  the  system  parameters  so  that  a proper 
weighting  can  be  given  in  the  filter  to  the  measurements. 
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The  form  of  the  modified  filter  is  thus  given  by 


the  equation 


where 


and 


£.  = x.  + G.  [z  - H.x  ] + /*. 
r i i i ii  i 


(3.11) 


x.  = <t>.  ..  .x.  + u. 

i i/i-l  i-l  i/i-l  i-l 


(3.12) 


X-  - filter  estimate  at  time  i after 
measurements 

x.  - predicted  state  at  time  i after 

Z . , measurements 
i-l 

- measurement  vector  at  time  i 

u^  ^ - input  vector  at  time  i-l 

G.  - filter  gain  modified  to  reflect 
parameter  uncertainty 

- filter  bias 


*i 


i'^i/i-l'**' i/i-l 


- matrix  relationships 
evaluated  with  or 


aQ  - expected  value  of  parameters  a 
The  bias  term  and  the  gain  in  the  filter  equation 
are  derived  in  the  following  subsections.  The  covariance 
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of  the  filter  error  is  also  determined.  In  addition,  the 
complete  set  of  algorithms  is  presented  and  a limited 
discussion  of  the  filter  stability  is  included  in  the  next 
section. 

3.2.1  Filter  Bias 

For  the  filter  to  be  unbiased,  the  expected  value 
of  the  difference  between  the  state  and  the  filter  estimate 
e^  must  be  zero  for  all  time  (i.e.,  Efe^]  = 0 for  i = 1, 

2,  3,...).  Utilizing  Equations  (3.4),  (3.10),  and  (3.11), 
this  difference  is  written 

a * 

e.  = x.  - x. 

ill 

= |"l  - G.H."U.  ..  .e.  _ 

L i ij  l/i-l  i-l 


+ [*  - - Vi  - "i 

+eL  k/i-i  - vA/i-ik-i 

k k 

+ JV-i  - Wi/i-ik-i 


+ a [^i/i-l  - 


5av  (3.13) 


where 


<\  ■ GiHi*i/i-l] 


a/i/i-1  - ViVi-1  - (3'14) 


[*i/i-l  - GiHA/i-l] 


- 'll...  , - G.  H , - G.H.  i b...  , (3.15) 

a.  i/i-l  la.  1 i/i-l  l iaji/i-l 


«k [ri/i-l  * 


/Vi-1  * Gi«i'Vi-1  - GlHia,-T'i/i-l  <3-16> 


and  where  for  simplicity,  (aQ) , («0^  and 

r.  , (a  ) have  been  shortened  to  f . </>...  , , and 
i/i-l  o l/l-l  l/l-l 

r • respectively,  and  the  second  order  terms  of  8a  have 
been  dropped.  Taking  the  expected  value  of  Equation  (3.13) 
and  recognizing  that  the  noise  inputs  are  zero-mean 
processes  and  are  uncorrelated  with  the  parameter  variations 
gives 
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e K]  - <i  - GiHi>  Vi-iECei-i] 


+ 11  - GiHi'ri/i.iEf“i-ii  - GiE[vi]  - ^ 


+ H " GiHi'*'i/i-l)  E[xi-l8“]J 


+ oj+i/ i-1  - GiHi',’i/i-l)  E[ui-l8o,J 


+ - GlHi/i/i-l)  E[wi-l8“k] 


-*i  +La<*i/i-l  - GiHiVi-l)  Etxi-l8“lJ 

k * 

(3.17) 


If  the  term  E[xi_i5ax]  is  approximated  by 


(3.18) 


The  bias  term  of  the  filter  is  then  given  by 


(3.19) 
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The  filter  error  can  thus  be  reduced  to 


ei  - (i  - GiV  Vi-i'i-i + (I  - GiHi>  - Givi 


k 


a.^i/i-l  ~ GiHi^i/i-l)xi-l 


+ „ (4>.  ..  . - G.H.  </».,.  ,)u.  , 
a,  l/l-l  1 1 l/l-l  i-l 


+ “v^i/i-l  ~ GiHi^i/i-l^  wi-] 


5 a, 


(3.20) 


This  expression  can  now  be  used  to  define  the  filter  gain. 
3.2.2  Filter  Gain 

Thfe  Kalman  filter  is  a minimum  error  variance  filter 
(see  Appendix  B) . It  is  therefore  desirable  to  achieve  this 
property  in  a modified  Kalman  filter.  If  the  gain  in  the 
modified  filter  is  selected  so  as  to  minimize  the  expected 
value  of  a quadratic  function  of  the  filter,  this  objective 
can  be  achieved.  The  equation  which  states  this  objective  is 


or,  alternately. 


(3.21) 
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1 I 


dtr  E[eieT] 


dG: 


= 0 


(3.22) 


where  tr  is  the  matrix  trace. 

Pursuing  the  above  objective,  the  covariance 
equation  must  first  be  derived.  Using  Equation  (3.20) 
and  recognizing  that  the  parameter  variations  are  uncor- 
related (i.e.,  EfSa^Sa^]  = 0),  the  covariance  may  be 
written 


T ,T 


+ (I-OiHi)ri/i.1E[wi_lw*.Jrf/i.1(I-OiHi)'r  + G^O^G* 


Ti  _T 


+ (I-GiHi)«i/l.1E[ei_1v,T_l]rT/i_i(I.GiH.) 


+ (I-SiHi)ri/i.1E[wi.1e^1]^/._i(i..0.H., 


- (I-GiHi)*i/i.1S[.i.1vJ]oJ  - OiE[vie^.1]^/i.1(I-OiHi) 


T 1AT 


- (I"GiHi)ri/i-iECui-iviJGi  - GiECvi"LiJrI/i-i(I-GiHi> 


T 1 r-T 


(cont) 
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— . ••• - 


k 


+ ak(<^i/i-l-GiHA/i-l)EtXi-leI-l8ak^  ^l/i-l“GiHiVi-l)T 


+ a^  ^ i/i-l-GiHi^i/i-l)  E[ui-lei-l®°k  ] ^i/i-l"GiHi^i/i-l) T 


+ (<*>.,.  ,-G.H  .<t>.  . ) E [e  . ,wT  ,8a.  ] (r /.  ,-G.H.T  ,)' 

i/i-l  i i i/i-l  L i-l  i-l  kJav  l/i-l  i li/i-l 


+ a (^i/i-l~GiHi^i/i-l)Etwi-lei-lSak]  ^i/i-l-GiHiA/i-l> 


+ [r.  --G.H  .r.  ,)e[w.  ,xT  ,8a-]  (<#>.,.  ,-G.H. <f>.  ..  ,)' 

i/i-l  l l i/i-l'  L i-l  i-l  kJa  ' i/i-l  l i^i/i-l' 

iC 


+ a ^i/i-l-GiHi^i/i-l)TE[xi-lwi-l8ak]  (ri/i-l“GiHi/i/i-l)T 
k 


+ ai  ^i/i-l~GiHi*^i/i-l^ ECui-lwi-l^Qfk]  ^i/i-l~GiHi^i/i-l^ T 

k 


+ ( /i/i-l“GiHi7i/i-l)  EtWi-lWi-l8ak^a  (ri/i-l*GiHi/i/i-l) 


(cont) 
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T T 

« (ri/i.i-GiHi/l/i-1)E[wi_1wi.1«fflk]  (Zi/i-i-GiHi/l/i.!) 
ak 

GiE  KxLl8ak]«kWVi-l‘GiHi*iA-l)T 

ov  <*i/i-l-GiHi*i/i-l>  E [xi-lvi5“  JGi 

- GiE[vtUi_16ak]akWi/i.1-GiHi+i/i.1)T 

- „ <*i/i-l-GiHi*i/i-l>E[ui-lvi-l6“lJGi 
k 

— T 

- GiE[viwi_1Sak]ak  (r L/i_l'GiHiri/i-l) 

- a (/i/^1-GiHiri/i.1)E[wi_1v^ak]G? 
k 

k 

+ ^i/i-l-GiHiVi-l)Etui-lui-lSak]a  ^i/ i-l“GiHi^i/ i-l) ^ 
ak  k 

T 2 T 

+ « (/i/i.i-GiHi/i/i.1)B[wi.1wi.18ork]  (ri/i_1-GiHiri/i.1) 
k x 

Qk  * 

(cont) 
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F 


+ - Vi-1}  E [ui-lxi-l8akL  (^i/i-rGiHA/i-l)T 


+ (t.,.  ,-G.H.  <fi.  ,)E[x.  _wT  (r.  /•  ,-G.H.r,.  ,)T 

ak  i/i-l  1 i i/i-l  L l-l  i-l  k-l«k  w i/i-l  i i i/i-l 


+ a ^i/i-i-GiHiri/i-i>EK-ixi-i8ak^ai  (*vi-i“GiHiVi-i)1 

k k 


(r ,.  ,-G.H.r,.  ,)e[w.  ut  so^]  ,-G.H.r.,.  ,) 

ak  i/i-l  i i i/i-l  u i-l  i-l  kJa^  i/i-l  i i i/i-l 


(3.23) 


The  expected  operations  in  this  equation  must  now  be 
evaluated.  The  term  E[ei_^e£_3j  is  simply  the  filter  error 
covariance  at  time  i-l,  which  will  be  designated  . The 

noise  inputs  w^_^  and  v^,  as  previously  defined,  are 
• uncorrelated  with  covariances  Qj__i  and  R^,  respectively. 
Terms  involving  single  functions  of  the  noise  inputs  are 
zero  since  the  terms  are  uncorrelated  with  the  noise  inputs 
and  the  expected  values  of  the  noise  inputs  are  zero.  An 
example  of  this  result  is 
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e[w.  nxT  .]  = E[w.  ] E [x . 

L i-l  i-lJ  L i-lJ  L i-lJ 


= 0 • E[xi_1] 


= 0 


(3.24) 


Terms  involving  the  noise  at  time  i with  results  at  time 
i-l  are  also  zero  since  future  noise  is  uncorrelated  with 
past  results.  Remaining  terms  containing  the  parameter 
variation  5a^  are  set  to  zero.  Other  terms  containing  5a^ 
are  estimated  using  the  variance  cr  . This  evaluation  and 
the  others  above  are  summarized  in  the  following  equations. 


E Te . .e™  ] = P.  , 

L i-l  i-lJ  i-l 


E [wi-iwLi3  ■ Qi-1 


Etviv3  = Ri 


E[e.  wT  ] = E[w.  eT  ] = 0 
^ 1-1  1-1J  U 1-1  1-1J 


E[e.  .vT]  = E[v.eT  .]  = 0 
L i-l  iJ  L i i-lJ 


(3.25) 


(3.26) 


(3.27) 


(3.28) 


(3.29) 
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Etwi-lvi]  * E tviwi- J = 0 

E[ei-ixLis“k] = E[xi-ieI-is«k^ = 0 
ECei-lui-lSak]  = E [ei-lui-ls“k]  = 0 

E[ei-lwi-l6aJ  = EK-lei-lSaJ  = 0 
EK-lxi-lSak]  * E fri-lwi-lSaJ  = 0 
EOi-lui-l8aJ  = El>i-lWi-lSakJ  = 0 
Et“i-iui-i8ak:i  = 0 
• E tvixI-i5“k;i  = Efxi-ivi6ak]  = 0 
Etviui-l8aJ  = E[ui-lviSaJ  = 0 


E [v . wT  5a.  ] = E[w  .v'J'Sa  ] = 0 


where 


= « i^i/i-1  + <3*49> 


(H.  = H.  <//.  /.  . + H. 

ak  i i/i-l'  ak  iaky i/i-l 


(3.50) 


(h ,r/.  .)  =/vH.r/.  . + h.  r. y.  , o.sd 

av  i i/i-l  a i i/i-l  iav  i/i-l 


The  modified  gain  is  therefore 


Gi  " Vi-l'w  (HA/i-l,T  +^/i-lOi-llHiri/i-l>T+Zavri 


li  + E lHiavEi  + *kHi°ksi' 


(3.52) 


where 


r.  - cr  d).  ..  . (P.  . +x.  AT  . ) (H . <£.  .) 

ak  i k ak  i/i-l  l-l  l-l  i-l  «k  ii/i-l 


+ ak^i/i-lui-lui-lak(Hi’/i/i-l)  + ak/i/i-lQi_lak(Hi/i/i-l) 


(cont) 


, ,i r fa* *■  r*”**  ■ 


+ ak^i/i-l^i-lui-lak(Hi^i/i-l)T 


? T 

+ <3-53) 


„ s.  = cr.  <£.  ,.  . (P.  +x.  .xT  ,)  (H.  <f>.  - )T 

«k  i k ^i/i-l ' i-l  i-l  i-l' ak  i^i/i-l' 


+ Wi-lUi-lak(HiVl)  +ri/i-lQi-lak(Hi/l/i-l)T 

+ 0.  ..  -X  uT  (H . tp  )T  + iff  u xT  (H . <f>  )T 
i/i-l  i-l  i-l  “k  l i/i-l  i/i-l  i-l  i-l  atk  i^i/i-l 

(3.54) 


This  gain  equation  contains  terms  found  in  the  Kalman  gain 
equation  and,  in  fact,  is  identical  to  the  Kalman  gain  when 
the  parameter  variations  are  zero.  Most  of  the  additional 
terms  in  the  expression  are  similar  in  form  to  the  nominal 
terms  of  the  Kalman  gain.  The  exception  is  the  quadratic 
term  of  the  input  u^  ^ and  the  cross  products  terms  of  the 
state  and  input  u^_^.  These  terms  occur  as  a result 

of  the  uncertainty  of  the  input  matrix 
3.2.3  Covariance  of  the  Filter  Error 

The  covariance  of  the  filter  error  at  time  is 
defined  by 


V * 


Pi  = E [(xi-xi)  (xi-xi)T]  = Efe^] 


(3.55) 


Using  Equation  (3.47)  gives 


Pi  * ( Vi-1'G i.Vi/i-11  Pi-1 ( A/i-l’ 


+ G 


. R . G . + V*  t . 
ill  i-*  at  i 


(3.56) 


where 


akti 


^♦i/i-l^iHiVi-l)  ^Pi-l+xi-lxi-l^ak  ^i/i-l“GiHi^i/i-l^ 


+ aklVi-l-GiHiVi-l)ui-lui-lok(Vi-l-°iHi*Vi-l,T 


(cont) 
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X * 


. 


(*4 


(3.57) 


The  first  three  terms  in  this  equation  correspond  to  the 
error  covariance  of  the  Kalman  filter.  The  additional 
terms  reflect  the  increased  uncertainty  in  the  estimate 
due  to  parameter  uncertainties  in  the  system  model.  Some 
of  these  added  terms  might  be  envisioned  if  an  expansion 
of  the  Kalman  filter  error  covariance  is  considered. 

However,  since  the  system  model  of  the  Kalman  filter  assumes 
a deterministic  input  u^  transmitted  through  a known  input 
matrix  ' terms  involving  the  input  are  not  con- 

tained in  the  Kalman  filter  error  covariance  (i.e.,  there 
is  no  uncertainty  in  known  quantities)  and,  therefore, 
would  not  be  anticipated  in  the  expansion.  Such  terms  do 
exist  because  of  the  parameter  uncertainties  within  the 
input  matrix  ^ ^ (a)  . 

3.2.4  Filter  Algorithms 

The  filter  algorithms  for  state  estimation  are 
summarized  here.  The  first  step  for  time  i-1  to  time  i 
is  the  computation  of  the  gain 


G. 

i 


J.  (H. J.  + R. 

ill  i 


+ 


(3.58) 
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X / 


where 


J.  = P.  ,H . + Y'  r . 
i i/i-l  i Z— iav  i 

k K 


s.  = y'  h.  s. 

1 ^ak  lak  1 


(3.59) 


(3.60) 


with 


(3.61) 


a aT  . , » T 

r.  = o-  <j>.  ,(P.  ,+x.  x.  ) (H.0  ) 

a i k a/i/i-1  i-l  i-l  i-l  av  1 i/i-l 

+ */i/i-iui-iuL«<Vi/i-i,T  vrvi-iQi-i«k(HiriA-i)T 


+ a ^i/i-lxi-lui-la  ^Hi^i/i-l^  +a  ^i/i-lui-lxi-la 
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(3.53) 


S.  = or  \<f>.  . (P.  +X.  JcT  ) (H.  <f>  ) 

ak  i fc|  i/i-l  i-l  i-l  i-l  <*k  i i/i-l 


+ Vi-iViui-iak(Vi/i-i,T  + yi/i-iQi-ie<k(Hiri/i-i)T 


(oont) 


<t> . , . .X  u'1'  (H  . ) 1 + 


Vi-1  i-1  i-lav  i i/i-1 
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The  bias  term  is  zero 


*L  = ° 


(3.1? 


The  filter  estimate  is  then  determined  similarly  to  the 
Kalman  filter  expression  in  Equation  (2.27)  by  the 
relationship 


(I  - G.H. ) x + G.z. 

ii  ii 


(3.6: 


where 


x . 

i 


~ <t>  3^ 

9i/i-l  i-1 


(3.6: 


Finally,  the  error  covariance  matrix  is  computer  from  the 
equation 


,(I  - G.H.) 


+ GiRiGi  + 


Ti 


(3.6 


(I  - G.H. ) P • /• 


+ (r,.  , G.H.r,.  jq.  , ( r -g.h  r. ..  .) 

i/l-l  l i l/i-l  i-la^  i/i-l  i 1 i/i-l 


+ (4>...  ,-G.H  ,)&.  -UT  . (<J/  . .-G.H>.  )T 

i/i-l  i i i/i-l  i-l  i-la^  i/i-l  i i i/i-l 


+ ('I'.  ..  -G.H  > _)u.  3CT  (0  -G.H>  ) 

a.  i/i-l  l i i/i-l  i-l  i-l«k  i/i-l  i i i/i-l 


(3.57) 


Note  that  all  matrices  in  the  expressions  above  containing 

the  parameters  a are  defined  assuming  the  nominal  parameter 

values  a . 

o 
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The  advantage  expected  for  this  filter  over  the 
Kalman  filter  is  the  higher  accuracy,  and  over  the  extended 
Kalman  filter  the  advantage  is  the  reduced  computation 
burden.  The  modified  filter  requires  only  n equations  for 
state  estimation  and  n(n+l)/2  equations  for  computation 
of  the  error  covariance.  With  p parameters  in  the  system 
equations,  the  extended  Kalman  filter  requires  n+p  equations 
for  state  estimation  and  (n+p) (n+p+l)/2  equations  for  the 
error  covariance  computation. 

The  gain  in  the  modified  filter  includes  additional 
terms  relative  to  the  standard  Kalman  gain.  In  general, 
when  the  parameter  uncertainties  exist  in  the  state 
equation,  the  effect  of  these  terms  is  to  increase  the  gain 
since  the  measurement  z^  is  a more  accurate  assessment  of 
the  state  than  the  prediction  of  the  state  given  by 
Conversely,  when  the  parameter  uncertainties  exist  in  the 
measurement  equation,  the  gain  is  decreased  since  the 
predicted  state  x^  is  a more  accurate  appraisal  than  that 
derived  from  the  measurement  z^. 

3.3  FILTER  STABILITY 

In  this  subsection,  the  asymptotic  property  of  the 
modified  Kalman  filter  error  is  considered.  The  approach 
utilized  in  the  assessment  has  been  applied  to  the  Kalman 
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filter.  However,  very  little  has  been  proven  or  is  known 
about  this  property  for  the  Kalman  filter.  For  this  reason 
and  because  of  the  difficulty  in  developing  general 
stability  results,  the  treatment  here  is  limited. 

The  assessment  makes  use  of  a stability  theorem 
commonly  defined  as  the  second,  or  direct,  method  of 
Liapunov.  This  theorem  can  be  stated  in  the  following 
manner.  Consider  a difference  equation  given  by 


= A yi-i 


(3.66) 


where  A is  a constant  matrix,  and  let  L be  a positive 
definite  Liapunov  function  of  y so  that  the  difference 
AL  is 


AL  = L (y±)  - L(yi_1) 


(3.67) 


Now,  if  AL  is  negative  semidefinite  (definite) , then  the 
difference  equation  is  stable  (asymptotically  stable) ; that 
is,  y^  remains  less  than  some  small  value  as  time  approaches 
infinite.  Conversely,  if  AL  is  positive  definite,  the 
difference  equation  is  unstable. 
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The  filter  error  as  defined  in  Equation  (3.20) 


1 


consists  of  a homogeneous  part  followed  by  terms  involving 

the  noise  inputs  and  perturbations  of  the  state  and  input. 

Since  the  terms  involving  the  noise  are,  on  the  average, 

zero,  they  can  be  ignored.  In  addition,  if  the  state  and 

input  remain  finite,  these  perturbation  terms  can  be 

ignored.  The  remaining  term  is  the  homogeneous  part  which 

% 

can  be  the  difference  equation  used  in  the  evaluation  if 
the  modified  gain  matrix  G^,  the  measurement  matrix  H^, 
and  state  transition  matrix  are  assumed  to  be 

constant  matrices  G,  H,  and  <f>,  respectively.  There  are 
many  problems  in  which  the  state  transition  and  measurement 
matrices  are  constant,  and  this  is  not  a major  restriction. 
However,  the  implication  of  a constant  gain  matrix  is  that 
the  filter  error  covariance  is  constant  (see  Section  3.2.4). 
In  many  problems,  the  covariance  would  not  remain  constant 
and,  in  fact,  might  be  increasing  in  value  with  increasing 
time.  Therefore,  this  evaluation  is  limited  to  determining 
values  of  the  gain  or  other  conditions  which  indicate 
stability. 

Therefore,  let  the  difference  equation  to  be 
evaluated  be  given  by 
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(3.68) 


yL  * (I  - GH)<^yi_1 

and  define  the  Liapunov  testing  function  as 

L = y^y-i  (3.69) 


which  is  clearly  positive  definite.  The  derivative  of  the 
Liapunov  function  is 

\ 


AT  T T 

al  = y.y.  - y._1yi_1 


T 

- GH)  (I 


(3.70) 


,T  T 

If  the  term  <P  (I  - GH)  (I  - GH)<f>-I  is  negative  semidef inite, 
then  the  filter  estimate  is  Liapunov  stable.  Equationally , 
this  is 


<*>T (I  - GH)T(I  - GH)4>-I  < 0 (3.71) 

The  modified  gain  can  assume  a range  of  values. 
Conceptually,  the  lowest  value  is  0 when  the  measurements 
are  disregarded.  In  this  case,  the  filter  is  stable  if 
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(3.72) 


4>t  4>-i  < o 

For  the  scaler  case,  this  is  equivalent  to 

<#>2  - 1 < 0 (3.73) 


or 


-1  < <i><  1 


(3.74) 


At  the  other  extreme,  the  measurements  may  determine  the 
state  exactly  (i.e.,  R^=0  and  a H^=0  for  all  i and  Tc)  . In 
this  case,  the  gain  G can  be  defined  by  (see  Appendix  C) 

GH  = I (3.75) 

The  derivative  of  the  Liapunov  function  then  becomes 


AL 


(“I)yi-1 
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(3.76) 


i 


which  is  unequivocally  negative  definite,  and  thus  the 
filter  estimate  is  stable. 

With  these  equations,  the  basis  is  provided  for 
developing  limited  guidelines  for  stability  with  the 
modified  Kalman  filter.  In  general,  the  filter  is  stable 
if 


4>T( I - GH)T(I  - GH)4>-I  < 0 (3.71) 

Specifically,  the  filter  would  be  expected  to  be  stable  if 
the  state  is  asymptotically  decreasing,  which  is  represented 
by  the  expression 


fy  - 1 < 0 (3.72) 

Additionally,  the  filter  is  expected  to  be  stable  if  the 
measurements  are  extremely  accurate.  This  occurs  when  GH 
is  approximately  equal  to  I.  For  other  conditions,  the 
actual  error  may  increase  with  time.  However,  the  error 
covariance  would  also  be  expected  to  increase.  The  filter 
error  is  Gauss ianly  distributed  and,  therefore,  the 
probability  of  any  such  error  should  be  adequately 
determined  with  the  filter  error  covariance. 
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CHAPTER  4 


ADAPTIVE  ESTIMATION  FOR  SYSTEMS 
WITH  UNCERTAIN  PARAMETERS 

The  Kalman  filter  derived  in  Chapter  2 is  am 
optimum  filter  for  stochastic  systems  with  known  para- 
meters. In  addition,  the  innovation  sequence  generated  by 
the  Kalman  filter  (i.e.,  the  sequence  of  differences  be- 
tween the  actual  measurements  and  the  predicted  measure- 
ments) is  a zero-mean  Gaussian,  white-noise  process.  When 
the  parameters  in  the  system  are  not  known  exactly,  however, 
the  filter  is  sub-optimum,  the  innovation  sequence  is 
altered,  and  adaptive  action  is  required  to  correct  the 
filter  performance.  In  this  chapter,  this  philosophy  is 
applied  to  state  estimation  for  systems  with  uncertain 
parameters.  An  adaptive  maximum  likelihood  Bayesian 
estimator  of  the  uncertain  parameters  is  derived  that 
functions  in  conjunction  with  a Kalman  filter.  The  Kalman 
filter  algorithms  operate  with  the  current  estimate  of  the 
parameters  to  generate  the  innovations  and  other  data 
utilized  by  the  Bayesian  estimator  to  improve  the  parameter 
estimates.  With  the  new  values  for  the  parameters,  the 
Kalman  filter  algorithms  are  again  operated  to  regenerate 
the  state  estimates.  Section  4.1  presents  the  basis  for 
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the  parameter  estimation  and  derives  a nonlinear  equation 
for  the  parameter  estimate.  The  solution  of  this  nonlinear 
equation  is  then  developed  in  Section  4.2.  The  adaptive 
filter  algorithms  are  summarized  in  Section  4.3,  and  the 
convergence  of  the  estimate  is  evaluated  in  Section  4.4. 

4.1  PARAMETER  ESTIMATION 

I 

This  section  develops  a Bayesian  parameter  estimate 
from  «the  variational  terms  of  the  conditional  probability 
density  function  of  the  parameters,  given  the  measurements. 
The  estimate  is  equivalent  to  a maximum  a posteriori 
estimator  known  also  as  the  maximum  likelihood  Bayesian 
estimator.  The  expected  value  of  the  variational  terms  of 
the  conditional  density  function  is  also  presented. 

4.1.1  Maximum  a Posteriori  Estimate  from  -In  p (a  Izn) 

One  of  the  most  useful  Bayesian  estimators  is  the 

maximum  a posteriori  estimator  which  is  obtained  by 

maximizing,  with  respect  to  the  parameters,  the  conditional 

density  function  of  the  parameters  p(a|zN)  where  ZN  is  the 

vector  of  realized  values  of  N measurements  z..,.z  . 

■1*  N 

Rather  than  directly  maximizing  the  conditional  density 
function,  however,  it  is  more  convenient  to  maximize  the 
negative  natural  logarithm  of  the  conditional  density 
function.  This  is  a valid  alternative  since,  for  any 
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function  p,  p and  -In  p attain  their  maxima  amd  minima, 

respectively,  at  the  same  point. 

Applying  Bayes ' rule  to  p (a  I Z ) results  in 

N 

p(Z  |a)  p (a) 

ptolV 2 (*•!) 

P'V 

Taking  the  negative  natural  logarithm  gives 

-In  ptalZjj)  = -In  p(ZN|a)  - In  p(a)  + In  p(ZN)  (4.2) 

Since  the  density  function  p(Z^)  is  not  a function  of  a, 
this  term  in  the  above  equation  will  not  affect  the  minimi- 
zation of  -In  p(a|^)  with  respect  to  a and,  therefore, 
need  not  be  evaluated.  The  remaining  terms  in  the  equation 
originate  from  the  joint  (unconditional)  density  p(a,Zjj). 
Consequently,  a derivation  initiated  with  the  density 
p(a,ZN)  would  produce  the  same  results.  Nevertheless,  the 
density  functions  p(Zjj|«)  and  p(ar)  must  be  evaluated  for 
development  of  a parameter  estimate. 

The  conditional  density  pCZ^a)  can  be  defined  in 


terms  of  the  density  function  of  the  innovation  sequence 
utilizing  the  derived  density  theorem  in  Appendix  D. 


r 


The  density  function  of  the  innovation  sequence  was  shown 
in  Chapter  2 to  be  Gaussian  with  zero  mean  and  covariance 
V.  for  an  optimum  filtfsr  with  known  parameters.  Using  this 
basis,  the  conditional  density  p(ZN|a)  is  defined  by  the 
equation 

p(ZN|a)  = P(JV*-I'n) 


V(2-)Nrn  ivj 


exp 


N m 1 

_%2>ivi  i,. 

i 


(4.3) 


i=l 


where  y.  are  the  innovations  and  | V^|  is  the  determinant  of 

T 

Vi  which  is  equal  to  HiPi/i.jHi  + The  density  function 

of  the  parameters  is  also  defined  to  be  Gaussian  and  is 

distributed  about  the  nominal  a with  covariance  cr.  The 

o 

parameter  density  is 


P («)  = 


Vu^l 


exp 


-Js(a-ao)To--1(a_Q:o) 


(4.4) 


where  j cr  | is  the  determinant  of  cr. 
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The  variationaii.  terms  of  the  negative  logarithm  of 
p(a|zN),  denoted  J(a),  can  now  be  determined  with  the  above 
equations.  Other  terms  in  -In  p(aJzN)  are  constant  and, 
therefoia,  will  not  affect  the  determination  of  the  maximum 
a posteriori  estimate.  The  quantity  J(a)  is  given  by  the 
equation 


N 


J<°>  -■£ 


j»Tv71f.  + lnlv.l 

1 1 i 1 


T -1 

(a-a  ) cr  (a-a)  (4.5) 

o o 


The  derivative  of  J(a)  with  respect  to  the  parameters  a is 
equivalent  to  the  derivative  of  -In  p(a|zN)  and  is 


dj(a) 

da 


i=l 


T 
V.\T 
ivi 


-1 


+ tr 


i da 


+ 2 cr  (a-a) 
o 


(4.6) 


where  tr  is  the  matrix  trace.  Note  that  the  partial 
derivatives  are  based  on  the  relationships  that,  for  any 
general  vector  a and  matrix  A which  are  functions  of  the 
parameter  a, 
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where 


aaTA-1a  _ „ daT  -1  T gAJ^ 
aS  2 dec  A a + a dQ[ 


(4.7) 


= -A-1  ^ A-l 
da  Oa 


(4.8) 


din  | A 1 
da 


dln|A|  d 1 A 1 

dUI  da 


i eUl 

A da 


tr  (A-1  |^) 
da 


(4.9) 


For  a minimum  value  of  dJ/da  and,  thus,  the 
maximum  a posteriori  estimate  of  the  parameters,  this 
equation  is  set  equal  to  zero 
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p 


This  equation  becomes 


N N 

E [j  (or)]  = 2tr(Ir)  + £lnlVJ  + tr(V 
i i 

N 

= Nr  + p + £ lnlvj  (4.12) 

i 

where  I is  an  identify  matrix  of  order  r which  is  the 
r 

dimension  of  the  innovation  vector,  and  similarly  I is 

P 

an  identity  matrix  of  order  p which  is  the  dimension  of 
the  parameter  vector.  This  equation  provides  an  assess- 
ment of  the  expected  value  of  J(a)  and  might  be  compared 
to  the  actual  value  of  J(a)  as  a check  of  the  estimation 
process . 

4.2  PARAMETER  SOLUTION 

The  derivation  of  the  previous  section  lead  to  an 
equation,  the  solution  of  which  yields  an  improved  estimate 
of  the  parameters.  However,  this  equation  does  not  have 
a general  closed  form  solution,  so  it  is  necessary  to 
consider  an  iterative  solution.  In  this  section,  an 
iterative  technique  is  developed  for  the  parameter  solution 
that  is  based  on  a method  known  as  the  "method  of  scoring" 
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in  the  statistical  literature  (see  Rao  1965) , the  modified 
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Newton-Raphson  method  in  the  control  literature,  and  also 
as  the  Gauss-Newton  method  in  other  contexts.  The  deriva- 
tives that  are  input  to  this  iterative  solution  are  also 
developed  in  this  section. 

4.2.1  Method  of  Solution 

One  means  of  obtaining  a solution  to  the  parameter 
equation  is  to  employ  the  Newton-Raphson  iteration.  The 
philosophy  of  the  Newton-Raphson  procedure  is  to  approxi- 
mate J(ar)  as  a quadratic  in  a by  using  a Taylor  series 
expansion  to  second  order  so  that 


J (<*) 


S J (a*) 


+ 


d J ( a *) 
da* 


T 

(a -a*) 


+ *5  (a  -a* ) 


a2J(or*) 


I da*‘ 


( a-a *) 


(4.13) 


where  a is  the  maximum  a posteriori  estimate  and  o*  is  the 
current  value.  The  derivative  is  then  approximately  a 
linear  function  of 


For  the  maximum  a posteriori  estimate,  the  term  dJ/da  is 
set  to  zero  which  provides  the  iterative  solution 


- 

-1 

• * 

d2J(a*) 

dJ  (a*) 

. &**2 

da* 

(4.15) 


The  determined  value  for  a is  then  a first-order  correction 
to  the  estimate  a*.  To  converge  to  a solution,  this 
equation  is  applied  repeatedly  until  the  corrections 
become  negligibly  small.  (That  is,  a*  is  used  to  calculate 
a.  This  a becomes  a*  for  the  next  iteration,  and  so  on.) 
Unfortunately,  the  determination  of  an  analytical  expres- 
sion for  the  second  partial  derivative  of  J(a)  in  the 
solution  equation  is  extremely  complicated  and  requires 
substantial  computational  time.  To  alleviate  this  diffi- 
culty, an  approximation  is  made  which  simplifies  the 
computation  while  maintaining  accuracy  over  large  samples. 

The  equation  for  J(a)  originates  from  the  varia- 
tional terms  of  -In  p(ZNla)  and  -In  p(a).  The  value  of 
the  second  partial  derivative  of  -In  p(ZN|ar)  is  approxi- 
mated by  the  expected  value  of  the  product  of  the  first 

partial  derivatives  of  -In  p(Z  | a)  and  is  known  as  the 

N 

"method  of  scoring"  in  Rao  (1965) . 
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a2ln  pCZjjIa) 


din  ptZjjIa)  I (din  p(ZN|a) 


Rao  has  shown  that  the  error  committed  by  using  this 
approximation  is  of  order  1/N  for  large  samples. 

The  approximation  is  motivated  by  development  of 
derivatives  of  the  probability  density  function.  A basic 
property  of  the  density  function  is  that 


go  go 

If  p(ZN|a)  dZN  = 1 


(4.17 


-00-00 


Performing  the  differentiation  with  respect  to  a parameter 
a and  using  the  fact  that  p is  equal  to  exp (In  p)  gives 

/V 


00  00  00  00 r 

— f--/p(ZNl<*)  dZN  =f"f  In  p(ZN|a)  p(ZN|or)  dZ^ 

aak4o_*S0  —oo  — "oo 


= e[-§-  In  p(zja)l 


jdinmi 


Continuing  the  differentiation  gives 


a 

da,  dc 


*■  ilY- 


— In  p(Z  | a)  p(Z  |a)  dz 

N J “ 


I a)  dZ 


+ /?V0L§_  In  p(Z>)l[jL  In  p(Z  |a)lp(Z 

-oo-ao  ^ak  N J[9«k  N 


2 

J In  p(zja)  + E.-&-  In  p(ZN|a)  In  p(ZN|a) 

^ j [K 


(4.19) 


From  this  equation,  the  developed  relationship  is 


'eW§S7  ln  PlZNl“’l  - 


A.  In  p (Z  | a)  A-  In  p(Z  | a) 

da.  N da,  N 


(4.20) 


The  second  derivative  of  J(«)  utilizing  this  result  is. 


therefore. 


8 


3 F (a,a)  + 2(t~18 
da2  k 


(4.21) 


where 


F(a,a)  = E 


dJAa)\/dljia) 


(4.22) 


with 


dHa ) 
da 


i 25l/i  -1 

^ ^ Vi  Vi 


T -1  SV±  _!  / . 

*iVi  iVi  ^i  + tr(Vi 


-1  ^i 


(4.23) 


and  where  F(a,o)  is  the  Fisher  information  matrix  and  8 

} 

is  the  Kronecker  delta  function. 

To  implement  the  iteration,  it  is  necessary  to 


evaluate  the  information  matrix.  Using  the  above  equations, 
realizing  the  independence  of  successive  innovations  for  a 
linear  filter  with  the  true  parameters  and  recognizing 
that  the  first  and  third  moments  of  a Gaussian  variable 
are  zero,  allows  the  expression 


8l 


w 


P(ak.«/)  * E 


* 


The  second  term  of  the  information  matrix  is  evaluated  in 


Appendix  E.  Utilizing  the  sample  data,  the  final  relation- 
ship is  then 


<W  ' % 


dv\  -1  dVi  ( 

4 i v.  — + 2tr(V 

da . 1 da . \ 


-i  av.  _!  avi 

* da  5a. 

k l 


(4.25) 


Equation  (4.25)  together  with  Equations  (4.21)  and  (4.15) 
thus  provide  the  parameter  estimate. 

4.2.2  Derivatives  of  Parameter  Solution 

The  parameter  solution  requires  the  input  of  two 
derivatives,  the  derivative  of  the  innovation  with  respect 


to  the  parameters  and  the  derivative  of  the  innovation 
variance  with  respect  to  the  parameters.  The  innovation 
is  given  by  the  equation 


The  derivative  of  is  readily  determined  for  any  defined 
problem.  The  derivative  of  the  expected  value  of  x^  can  be 
developed  from  the  equation 


xi = Vi-A-i + *i/i-iui- 


i/i- 


^i/i-1  xi-l  + Ki-lVi-l  + ^i/i-lui-l  (4 . 27) 


which  gives  the  result 


dx^^  d(f> 


i/i-l  ft  + u 

da  da  da  i~i 


+ Vi-I 


+ 1 + K.  ^ 

da  da  da 


(4.28) 


This  equation  demonstrates  a recursive  behavior  in  the 
development.  To  complete  and  evaluate  the  equation,  the 
derivative  of  the  gain  must  be  defined.  This  derivative 
is  determined  from  the  equation 


% X 


1 


Ki-1  - Pi-l/i-2HI-l  (Hi-lPi-l/i-2HI-l  + Ri-Irl 


* W^I-i 


(4.29) 


and  is  given  by 


**i.l  9Pt-X/i-2llT  V-1  +p,  ?fkvTl 

da  da  i-1  i-l/i-2  da  i-l 


- Pi-l/i-2Hi-lVI-l 


SVj-1 

da 


viii 


•;  /i  o T -1 

—J-?  1-~-  H ,V.  . + P.  ...  _ 

3a  l-l  i-l  1-1/ 1-2 


ffk  v-1 

da  i-1 


aVi-l  -1 
- Ki-l  Vi., 


(4.30) 


The  derivatives  of  the  conditional  state  error  covariance 

N ■( 

and  the  innovation  variance  must  now  also  be  determined. 

The  derivative  of  the  conditional  state  error 
covariance  is  developed  from  Equation  (2.32)  to  yield 
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— - ....  j 


a^i_-l/i-2  _ a<^L-2/i-2  T 5Pi-2  fT 

da  da  *-2*1-1/ i-2  + ^i-l/i-2  da  ^i-l/i-2 


+ <f>  p (d<t>i-l/i-2^+  ri-l/i-2  0 rT 

®-l/i-2Pi-2  1 da  / 5a  Qi-2/i-l/i-2 


+ ^"i-1/  i-2^i-2 


'i-l/i-2 


(4.31) 


The  derivative  of  the  state  error  covariance  is  now 
required.  Using  Equation  (2.26),  this  derivative  is 


aPi-2  ‘aKi-2 


da  Hi-2Pi-2/i-3  " *i-2  da  *i-2/i-3 


aiii-2 


+ (I  - K H ) aPi~2/i-3 
(I  Ki-2Hi-2J  d a 


(4.32) 


(If  the  parameter  occurs  only  in  the  input  matrix  «/»,  this 
derivative,  the  conditional  error  covariance  derivative, 
and  the  gain  derivative  are  zero.)  These  expressions 
produce  a recursive  relationship  for  the  derivative  of  the 
innovation.  The  relationship  requires  derivatives  from  the 
penultimate  time.  Computationally,  however,  the  numerical 
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procedure  cam  be  developed  so  that  only  derivatives  from 
the  preceding  time  need  to  be  saved. 

The  other  derivative  needed  in  the  parameter 
solution  is  the  innovation  variance  derivative.  This 
derivative  is  developed  from  the  equation 


Vi  = HiPi/i-lHI  + Ri 


(2.44 


The  derivative  is 

av . aH . ap . / . , aHT 

-dt  = -d  Pi/i-lHi  + Hi  a,1'  Hi  + HiPi/i-l  ST  (4-33 

and  is  zero  when  the  parameter  occurs  only  in  the  input 
matrix  The  derivative  of  in  the  equation  is  readily 
determined  while  the  derivative  of  is  defined  by 

the  recursive  relationship  of  Equations  (4.30),  (4.31), 
and  (4.32) 

4.3  SUMMARY  OF  ADAPTIVE  FILTER  ALGORITHMS 

Previous  sections  provided  the  basis  for  improving 
the  parameter  estimates  and,  therefore,  allowing  a more 
accurate  state  estimate.  The  recursive  equations  of  the 
filter  are  summarized  in  this  section.  The  initial 
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A / 


algorithms  parallel  the  Kalman  filter  in  computation  and 


require  storage  of  the  time  histories  of  matrices  ^(a), 

-fi/i.i  («)  * H^a),  Qi#  and  and  of  the  computed 

matrices  at  each  step  x.  , x.  , K. , P.  and  P.  ..  , . The 

i i i i i/i-l 

initial  filter  algorithms  at  time  are 


xi  ■ Vi-A-i  + 


(4.34) 


P.7.  » <t>.  ..  .P.  >T/.  . + F /•  .Q.  .A  . (4.35) 

i/i-l  i/i-l  i-l  i/i-l  i/i-l  i-l  i-l 


Vi  = HiPi/i-lHi  + Ri 


Ki  ” Pi/i-lHiVi1 


x.  = (I  - K.H. )x.  + K.z. 

l ill  ii 


(2.44) 


(4.36) 


(2.27) 


T T 

P.  = (I  - K . H , ) P . ..  _ (I  - K.H.)  +K.R.K.  (2.28) 

i ii  i/i-l  il  ill 


The  initial  conditions  for  the  computation  are  noted  here 


as 


xi=o  = xo 


(4.37) 
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I * 


p . = p 

o/ — 1 o 


(4.38) 


V = 0 
o 


(4.39) 


Kc  = 0 


(4.40) 


Rc  = 0 


(4.41) 


Note  also  that  Q is  not  zero,  but  is  defined  by  the 

o 


numerical  problem. 

To  enable  an  estimate  of  the  parameters,  deriva- 
tives of  the  Kalman  filter  terms  must  also  be  computed  at 
each  time.  The  first  required  derivatives  are  those  of 
the  conditional  state  estimate  and  the  conditional  error 
covariance  which  are  determined  prior  to  the  Kalman  filter 
algorithms  by  the  equations 


da 


ak  xi-l  da  u 


1-1 


+ <t> 


i/i-1 


dxi-l 


dK 


i-1 


da 


da 


"i-l  + Ki-1 


d"i-l 


da 


(4.28) 
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aPi/t-l  «*i/i-l  .T  j.  aPi-l  aT 

gZ ~~dt Pi-A/i-l  + *1/1-1  da  *1/1- 


64?  , dT.  ..  ..  . t 

+ — £r+  al1’-  Qi/lfi/i-1 


+ /7  /;  i Q.; 


a^/i-l 


i/i- 

i/i-lwi-l  da 


The  initial  conditions  for  these  equations  are 


(4.42) 


(4.43) 


(4.44) 


(4.45) 


(4.46) 


The  next  derivatives  are  computed  after  the  Kalman  filter 
algorithms.  These  derivatives,  in  the  order  of  their 
computational  sequence,  are 


VHVPBVNMPSP'IOTPI 


dx. 


dv i eHt  _ _ 

da  i i 5a 


5a 


(4.26) 


5Vi 


5H. 


5T  pi/i-lHi  + Hi  -0^~  Hi  + HiPi/i-l  5^  (4'33) 


5k.  5p.  /.  , m 

1 iA-1 


5H 


5a 


aa  Hivi  + pi/i-i  ~dt  vr  - Ki  a^r  vr  <4-47> 


i ..-1 


Vi  - K, 


avi  -1 
»»  * 


api 


da 


a k aH. 

HiPi/i-l  " Ki  pi/i-l  + (J  - KiHi^ 


ap 


i/i-l 


5a 

(4.48) 


The  accumulated  score  S . and  the  information 

i 


Pi^ak'a/)  are  also  determined  at  each  time  point  by  the 
relationships 


S.  = 


„ T 

dvi  1 
2 -3-i-  V.J/. 
da  i i 


5V, 


- V 


i i da  i 


T vvi  i dv i 

tV.1  V.1*/.  + tr  (V . 1 -^-i)  + S. 
' i 5a 


i-1 


(4.49) 
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Fi  <“*•“/> 


Ok?  -i  du,  / .1  Sv  x 6v  \ 

- 4 &r  v i -sr  + 2 tr(vi  tst  vi  -sr) 


+ Fi-i(V“/ 


(4.50) 


These  quantities  are  initially  zero  and  are  determined 
only  every  N samples  so  that,  when  i is  equal  to  N,  an 
estimate  of  the  parameters  is  made  and  the  quantities  are 
reset  to  zero.  The  parameter  estimate  a is  made  by  the 
relationship 


a = a*  - 


a2J(a*) 

- i 

dJ  (or* ) 

. da*2  . 

. da*  . 

(4.15) 


where 


g-..J. (■?!)-  = fm  (a  * , a* ) + 2o7  \ 


da*2 


N'“  7 ’ ""'k  6k f 


(4.51) 


2^1  . S„  + 2<r-h«*-0) 


(4.52) 


and  where  or*  is  the  current  parameter  estimate.  After  the 
estimate,  the  last  N measurements  are  reprocessed  with  the 
algorithms  to  establish  a new  parameter  estimate.  This 
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process  is  repeated  until  the  corrections  to  the  parameter 
estimate  become  negligibly  small.  At  this  point,  the 
state  estimates  produced  by  the  filter  are  accepted. 

Also,  after  each  cycle  of  N measurements,  derivatives  are 
reinitialized  so  that 


(4.53) 


(4.54) 


(4.55) 


= 0 


(4.56) 


This  procedure  is  termed  batch  processing.  The  data  from 
the  previous  batch  is  considered  as  a priori  information 
for  the  current  batch  and  so  on  until  the  estimates  are 
accepted. 

4.4  CONVERGENCE  OF  ESTIMATE 

Consider  now  the  convergence  of  the  parameter 


estimate.  The  estimate  is  defined  with  the  derivative  of 
J(a)  set  to  zero  or 


dJ(q) 

da 


= 0 


(4.57) 


This  derivative  can  be  redefined  as 


N 


d>J  (a)  _ , _ , _ , 

~dZ~  * aJ(Va) 


(4.58) 


Now  let  the  expected  value  ffj(z  ,a|a  ) of  the  individual 
terms  aJi(zi»a)*  conditioned  on  the  true  parameters,  be 
defined  as 


•00  00 


=J  JaJi(zi.a)£(zi\at)  dz±  (4.59) 


-00  -00 


where  a if  the  true  value  of  the  parameter,  so  that  the 
average  of  the  individual  values  if 


N 


„J<ZK.a|at) 


(4.60) 
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Therefore,  aJ(ZN<a)  is  the  sample  mean  from  a population 
whose  mean  is  a J(ZN,a|  «t)  . By  the  weak  law  of  large 
numbers,  converges  in  probability  to  aJ(ZN,a|at) 

which,  for  an  arbitrary  value  of  € greater  than  zero, 
is  written 


lim  P 
N — ►<» 


aJ(ZN,a)  -aJ(ZN,a|at)  > c 


= 0 


(4.61) 


Since  the  parameter  estimate  is  based  on  minimizing  J(o), 
the  mean  value  aJ  (Z^,o:  | °^)  is  strictly  increasing  over  a 
range  of  a including  a^.  „Let  the  range  be  - 8 to  a + 
where  5 is  some  positive  value.  If  a is  set  equal  to  afc, 
then 


5 


„J(Z  , a la)  = E 
« N t'  t 


= 0 


(4.62) 


Consequently,  when  a is  a -8,  the  value  of  J(Z„,<*|qJ  is 

t a N 1 t 

given  by  the  inequality 


aJ<W*lat>  * 0 


(4.63) 
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CHAPTER  5 

X If 


NUMERICAL  RESULTS 

This  chapter  applies  the  filter  estimation  techni- 
ques of  the  previous  chapters  to  a numerical  example.  The 
example  consists  of  estimating  the  attitude  rates  of  a 
symmetric  spinning  body,  assuming  an  inexact  knowledge  of 
the  spin  rate.  Results  are  derived  for  both  the  modified 
Kalman  filter  and  the  adaptive  filter  presented  in  the 
previous  chapters.  Numerical  results  are  also  generated 
for  the  Kalman  and  extended  Kalman  filters  to  provide  a 
comparative  assessment.  The  numerical  example  and  the  four 
filter  formulations  are  delineated  in  the  first  sections  of 
the  chapter.  The  last  sections  provide  the  numerical  re- 
sults and  comparisons. 

5.1  NUMERICAL  EXAMPLE  , 

The  angular  motion  of  a moment-free  spinning  body 
is  described  by  a set  of  equations  known  as  Euler's  equa- 
tions. These  rotational  equations  of  motion  can  be  derived 
by  applying  the  principle  of  conservation  of  angular 
momentum.  That  is,  the  time  derivative  in  inertial  space 
of  the  angular  momentum  vector  determined  about  the  system 
mass  center  is  equal  to  zero.  For  a symmetric  body 
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spinning  at  a constant  rate  about  the  principal  axis  of 
maximum  moment  of  inertia,  designated  here  as  axis  3,  the 


moment-free  equations  of  motion  are: 


(5.1) 


(5.2) 


(5.3) 


and 

- angular  rate  about  body  axis  1 of  an 
orthogonal  set  of  body-fixed  axes 

w2  - angular  rate  about  body  axis  2 of  an 
orthogonal  set  of  body-fixed  axes 

w3  - angular  spin  rate  of  body  (assumed  to 
be  constant) 

1^  - principal  moment  of  inertia  about  axis  1 
(also  moment  of  inertia  about  body  axis  2 
for  symmetric  body) 

1 2 - maximum  moment  of  inertia  for  principal  axis  3 

A - moment-of-inertia  factor 


ti>^  + Acti2a*3  = 0 


a»2  - A&>^a»3  - 0 


where 


A = —2.  - 1 

I, 
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These  equations  can  be  defined  by  a discrete  representation 
utilizing  Appendix  A.  Assuming  a short  sample  period 
relative  to  the  natural  frequency  of  motion  and  allowing 
disturbances  in  the  form  of  additive  white  noise,  the 
equations  of  motion  become 


V 

“2. 

i+1  *■ 

L<u3i 


3 

1 


"l' 

+ 

W1 

J 

o>2 

i 

.W2. 

(5.4) 


where 


r = XT 


(5.5) 


and  where  T is  the  sample  period  and  w is  the  disturbance 
noise. 

Spin-stabilized  spacecraft  typically  spin  at  60  rpm 
or  6.3  rad/sec  and,  because  of  volume  and  packaging  con- 
straints, have  moment-of-inertia  factors  X of  approxi-  , 
mately  0.1  [e.g.,  Williamson  (1974)  reports  a spin  rate  of 
60  rpm  for  the  FLTSATCOM  satellite  and  a moment-of-inertia 
factor  of  about  Q.l].  Therefore,  let  the  nominal  value  for 
^ be  6.3  rad/sec  and  let  X be  0.1.  Also,  let  the  sample 


OJ 
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period  be  0.1  sec  so  that  t is  0.01  sec  and  o^t  is 
nominally  0.063.  Let  the  variance  of  the  spin  rate  be 
0.3  (rad/sec)  and  set  the  true  value  of  the  spin  rate  at 
7.4  rad/sec  (i.e.,  a value  roughly  equivalent  to  twice  the 
standard  deviation) . Consequently,  the  true  system 
dynamics  are  given  by 


OU-i 


LW2J 


i+1 


1 - 0.074 

0.074  1 


"1 


w. 


LW2J 


(5.6) 


and  the  assumed  system  dynamics  are 


1 - 0.063' 

0.063  1 


“1 

Ci)  2 

i+l 

"l 

oj2 


Ji 


w. 


Lw2 


(5.7) 


Let  the  statistics  on  the  initial  conditions  be 
defined  as 


r a 

” * 

i ^ 

0.63 

CM 

<3 

n 

. 0 - 

rad/sec 


(5.8) 


"NT- 


M0  ■ 


O.OOOl  0 

0 0.0001 


(rad/sec) 


(5.9) 


and  let  the  actual  initial  conditions  be  set  at  approxi- 
mately 3cr  values  in  order  to  observe  the  transient  behavior 
of  the  system,  or 


* 

1 

0. 66*1 

1 

CM 

3 

— i 

r\ 

o 

• 

o 

u> 

I 

rad/sec 


(5.10) 


Assume  the  noise  covariance  of  the  system  to  be  time 
invariant  and  given  by 


Qi  = 0.00001 


1 0 
0 1 


(rad/sec) 


(5.11) 


Let  the  measurements  be  related  to  the  angular 
rates  with  additive  white  noise  by  the  equation 


‘2J 


1 0 
0 1 


LW2J 


lV2j 


(5.12) 
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I 


and  assume  the  covariance  of  the  measurement  noise  to  be 
time  invariant  and  given  by 


Ri 


0.0005 


1 0 
0 1 


2 

(rad/sec) 


(5.13) 


Thus,  these  relationships  now  define  a numerical  example 
that  can  illustrate  the  filter  estimation  techniques  of 
the  previous  chapters.  In  addition,  the  assumed  values 
are  consistent  with  physical  constraints  of  spinning  bodies 
identified  in  Appendix  F. 

5.2  KALMAN  FILTER  FORMULATION 

The  state  to  be  estimated  by  the  Kalman  filter  is 
designated  and  given  by 


(5.14) 


The  predicted  value  of  the  state  based  on  the  nominal  values 
of  the  transition  matrix  is 


A 
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1 

0.063 


0.063 

1 


rA 

x. 


1*2 


i-1 


where  the  error  covariance  is  given  by 


P . . . , = 4>.  .P.  10T/.  . + Q . . 

l/l-l  1/1-1  1-1 'l/l-l  1-1 


1 

- 0.063 

1 

- 0.063" 

— 

0.063 

1 . 

Pi-1 

0.063 

1 

0.00001  0 

0 0.00001 

m 


The  initial  values  input  to  filter  predictions  are 


po  = 


0.0001  0 


0 0.0001 


(5.18) 


The  gain  equation  for  the  filter  estimate  is  given  by 


Ki  - pi/i-i  [pi/i-i  + Ri' 

m 


-1 


= P. 


i/i-1 


P. 


i/i-1 


0.0005  0 


0 0 . 0005/ 


1-1 


(5.19) 


The  filter  estimate  is 


, = (l-K.lx.  + K.z. 

i L ij  i i r 


(5.20) 


where  the  measurements  z ^ are  derived  from 


LZ2 


Lx2 


LV2J 


(5.21) 


and  the  state  is  determined  from  the  true  system 
dynamics 
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1 - 0 . 0741 


(5.22) 


0.074 


with  the  true  initial  conditions 


X1 

0.66 

-X2  - 

r\ 

0.03, 

0 


(5.23) 


Finally,  the  covariance  of  the  filter  error  is 


P.  = fl-K.lp...  fl-K,  ]T  + K.R.KT 
i L ij  i/i-l L ij  ill 


r l r it 

0.0005 

0 

LI~KiJpi/i-iLI-KiJ  +Ki 

. 0 

0.0005 

a 

(5.24) 


These  quantities  determined  at  time  i become  the  inputs  to 
the  filter  equations  for  the  estimate  of  the  state  at  time 
i+1.  Thus,  the  equations  form  a recursive  set  of  expres- 
sions that  process  the  measurements  as  they  become  avail- 
able to  the  filter. 
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5.3 


MOD IFIED  KALMAN  FILTER  FORMULATION 


The  modified  Kalman  filter  operates  in  a similar 
manner  to  the  Kalman  filter  with  the  exception  that  the 
gain  and  error  covariance  matrices  are  altered  to  reflect 
the  uncertainty  in  the  system  parameters.  The  modification 
in  these  matrices  is  accomplished  by  utilizing  the  para- 
meter covariance  and  the  derivative  of  the  system  matrices 
with  respect  to  the  parameter.  For  this  example,  the 
uncertain  parameter  is  the  body  spin  rate  a^.  The 
derivative  of  the  transition  matrix  with  respect  to  the 
spin  rate  is 


a^i/i-1 


0 -0.01' 
0.01  0 


(5.25) 


2 

and  the  spin  rate  covariance  cr  is  0.3  (rad/sec)  . The 
terms  of  the  gain  matrix  become 


<x^ri  °a^i/i-l  Pi-1  + Xi-lXi-lj  a^j^i/i-lj 


0 

- 0.01 

^ A AT  1 

pi-l  + xi-lxi-l 

0 

- O.Ol' 

0.01 

0 

0.01 

0 

(5.26) 
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si  = 0 


(5.27) 


p=l 

Ji  = Pi/i-l  + X “kri 

k=l 


" Pi/i-l  + °*3 


0 - 0.01 
0.01  0 


[pi-l + 


0 - 0.0l'|T 

0.01  0 

(5.28) 


and  the  gain  matrix  is 


G.  = J. (J.  + R.  + S. ) 

i ii  i i 


-1 


Ji 


Ji  + 


'0.0005  0 

0 0.00057 


-l-l 


(5.29) 


The  filter  estimate  is 


Xi  = (I-Gi)xi  + Gizi 


(5.30) 


and  the  error  covariance  of  the  filter  estimate  is  given  by 


pi  - (I-Gi>pi/i-l(I-Gi)T  + GiRiGI  + Ti 


0.0005  0 

= (I-Gi)pi/i-l(I-Gi)T  + Gi  Gi  + Ti  (5-3D 

0 0.0005 


where 


P=1 

■ £ «v‘i 


A AT 


0.3  (I-Gi) 


0 - 0.01 


0.01 


{Pi-l  + ^i-l^i-l) 


0 - 0.01 


0.01 


(I-G±) 


(5.32) 


Note  that,  in  this  filter,  the  state  estimate  uses  one  less 
equation  and  the  error  covariance  matrices  require  five 
less  equations  than  the  extended  Kalman  filter. 
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X 0 


5.4 


EXTENDED  KALMAN  FILTER  FORMULATION 


In  the  extended  Kalman  filter,  the  uncertain  para- 
meter in  the  system  dynamics  is  augmented  to  the  state 
vector  so  that 


The  system  dynamics  are  then 


(5.33) 


(5.34) 


From  the  development  in  Chapter  2 , the  predicted  value  of 
the  augmented  state  is  written 


x£  - 0.0 lx^ 

x^  + 0. Olx^x^ 


(5.35) 


and  the  covariance  of  the  predicted  state  becomes 


i 


1 - 0.01x1  - 0.01x1 

3 2 


= 0.01x1 


0 . 01x£  P|  . 


1 - 0.01x1  - 0.01x1 


0.01x1 


0 . Olx^ 


0.00001 


+ 0 0.00001  0 


(5.36) 


The  initial  estimates  provided  to  the  filter  are 


1 0 . 63 i 


A, 

Xi  = 


(5.37) 


p; 


0.0001  0 0 


0 0.0001  0 


0 0.3 


(5.38) 


The  gain  equation  of  the  extended  filter  is 


Ki  - pi/i-iHi(Hipi/i-iHi  + V 


-1 


= P!  „ 


i/i-1 


10  0 
0 10 


10  0 
0 10 


pi/i-l 


10  0 
0 10 


T 

+ 


0.0005  0 

0 0.0005 


-1 


(5.39) 


The  filter  estimate  using  this  gain  is 


» x'  + K' (z  -x’ ) 

i i ill 


(5.40) 


and  the  error  covariance  of  the  estimate  is  given  by 


T T 

p:  = (i-k:h:)p'  . ,(i-k:h:)  + k:r.k: 

i 11  i/i-l  11  ill 


i-k: 

i 


1 0 0] 
0 10 


pi/i-l 


I"Ki 


10  0 
0 10 


+ K! 
1 


0.0005  0 


Ki 


0 0.0005 < 

(5.41) 
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The  measurements  are  given  by  the  equation 


V 

1 0 o' 

r 1T 

vi" 

= 

|xix2x3j  + 

-Z2- 

. 

0 10 

L -J 

-V2- 

(5.42) 


where  the  state  is  based  on  the  true  system  dynamics 


x ' 

1 

X' 

— 

2 

X* 

L 3. 

i+1 

- 0.074  0 
1 0 
0 1 


x ‘ 

w. 

1 

1 

X 1 

+ 

w„ 

2 

2 

7.4 

0 

i 

. 

i 

- 

(5.43) 


with  the  true  initial  conditions 


x’  ' 

'0.66' 

1 

X' 

— 

0.03 

K) 

x- 

7.4 

L 3. 

» 

(5.44) 


The  extended  Kalman  filter  operates  with  the  above  equations 
in  a recursive  manner  estimating  both  the  states  and  the 
.uncertain  parameter  of  the  system  dynamics. 
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5.5 


ADAPTIVE  FILTER  FORMULATION 


The  adaptive  filter  uses  the  algorithms  of  the 
Kalman  filter  and  derivatives  of  the  filter  terms  to 
perform  the  parameter  estimation  and  improve  the  state 
estimates.  The  filter  algorithms  operate  with  the  current 
estimate  of  the  body  spin  rate  designated  a*.  The  initial 
value  assumed  for  the  spin  rate  is  6.3  rad/sec.  The 
process  is  initiated  by  computing  the  two  derivatives 


a*. 


i/i-l  a 

da  Xi-1  + *i/i-l 


a*i-l 


dK 


i-1 


da 


da 


*i-l  + Ki-1 


dv  i_i 


da 


p 

- 

_ 

“ 

0 

- 0.01 

A 

1 

- 0.01a* 

— 

X.  . + 

0.01 

0 

1 — X 

O.Olot* 

1 

6x 


i-1 


dK 


i-1 


dv 


da 


da 


Vi-1  + Ki-1 


i-1 


da 


(5.45) 


api/i-i 

da  da 


’i/i-l  ,T 

Pi-l*i/i-l  + 


dP.  . 

• J--1  uT 

Vi-1  da  0i/i-l 


T 

d<t> 

— j/.i-l 

da 


(cont) 
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' 0 

- 0.0l' 

1 

- 0.01a*' 

= 

0.01 

0 

pi-l 

0.01a* 

1 

1 

- 0.01a*‘ 

3Fi.i 

’ 1 

- 0.01  *' 

0 . Ola* 

1 

da 

0.01a* 

1 

1 - 0.01a*" 

r 

rH 

O 

• 

o 

1 

o 

1 

+ 

P • i 

0 . 01a*  1 

i-l 

0.01  o 

(5.46) 


The  initial  values  of  the  derivatives  in  these  equations 
are  zero.  In  addition,  these  derivatives  are  set  to  zero 
after  processing  every  N measurements  for  the  parameter 
estimate.  In  this  example,  N is  set  equal  to  30. 

Following  computation  of  the  above  equations,  the 
Kalman  filter  algorithms  are  utilized.  These  algorithms 
are 


Vi-1* 


i/i-1 


1 

0.01a* 


(5.47) 


+ Qi-1 


1 - O.Olct*! 


0 . Ola* 


Pi-1 


1 - 0.01a* 

0 . Ola*  1 


0.00001  0 


0 0.00001 


Vi  “ Pi/i-l  + Ri 


= P.  . + 
i/i-1 


0.0005  0 

0 0.0005 


K. 

i 


P 


i/i-1 


-1 

V. 

1 


= p 


i/i-1 


'0.0005 


P.  ..  . + 
i/i-l 


0.0005. 


= (I-K.)x.  + K.z. 
^ ^ « « « 


L2 


= (I-Ki)Pi/i.1(I-Ki)T  + KiRiKT 


0.0005  0 

= (I-K.JP.-  I (I-K-) T + Kl  kT  (5.52) 

i i/i-i  i 0 0.0005 


where  the  initial  conditions  and  measurements  are  identical 
to  those  of  the  Kalman  filter  presented  in  a previous 


section. 


The  following  derivatives  are  also  computed  at  each 


time  sample 


dvi 


(5.53) 


aVi  aPi/i-l 


(5.54) 


dKi  _ apiAdL  v-i  . K v-1 

da  i i da  a- 


i/i-] 


-et  - - -5T  Pi/i-l  + da 


(5.55) 


(5.56) 
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These  derivatives  form  the  accumulated  score  and  information 


determined  by  the  equations 


dv]  , m i av.  , 

s.  = 2 vT  v.  + v v~  vT  i/. 

l 5a  ii  l l 5a  i i 


+ tr 


V 


-l  dvt 
da 


+ Si-1 


F.  = 


a T 

dv  _L 

4 —A  V.  — i + 2 tr 
5a  1 da 


ei,i 


-i  -i 

V.  =•  V.  i 

i da  1 da 


+ F 


i-1 


The  parameter  estimate  is  then  determined  by 


a = a*  - 

F + 2<r~1 

-1 

S„  + 2cr“1(a*-a  ) 

N 

N o 

= a*  - 


F + 

FN  + — 


-1 


SN  + ^ (a*-6. 3) 


(5.57) 


(5.58) 


(5.59) 


After  the  parameter  estimate,  the  last  N measurements  are 
reprocessed  repeatedly  until  the  correction  to  the  para- 
meter estimate  is  negligibly  small  (0.01  in  this  example) . 

At  this  point,  the  state  estimates  are  accepted  and  the 
next  series  of  measurements  is  processed.  To  summarize,  the 


Kalman  filter  algorithms  continue  to  estimate  the  two  sys- 
tem states  using  the  current  estimate  of  the  uncertain 
parameter  while  the  additional  algorithms  adaptively  deter- 
mine the  uncertain  parameter  from  the  innovations. 

5.6  NUMERICAL  COMPARISON 

The  state  estimates  and  the  errors  in  the  estimates 
are  graphically  displayed  in  the  following  figures. 

Figures  1 through  4 show  the  estimate  of  state  1 compared 
with  the  true  value  of  the  state  for  the  Kalman  filter,  the 
modified  Kalman  filter,  the  extended  Kalman  filter,  and  the 
adaptive  filter,  respectively.  Both  the  estimate  and  the 
true  value  of  the  state  have  been  normalized  to  the  maximum 
value  of  the  state.  In  Figures  5 through  8,  the  percentage 
difference  between  the  normalized  estimate  and  the  state 
for  each  filter  is  presented.  The  Kalman  filter  produced 
a maximum  error  in  the  estimate  of  about  8 percent  while 
the  modified  Kalman  filter  reduced  the  maximum  error  to 
about  6 percent.  By  estimating  the  uncertain  parameter  in 
the  system  dynamics  (i.e.,  state  3),  the  extended  Kalman 
and  adaptive  filters  dramatically  reduced  the  error.  The 
maximum  error  of  the  extended  Kalman  filter  was  approxi- 
mately 2 percent  while  the  adaptive  filter  was  about  the 


same  but,  at  certain  times,  slightly  larger.  Similar 
filtering  results  for  state  2 are  presented  in  Figures  9 
through  16. 

Only  the  extended  Kalman  and  adaptive  filters 
provided  estimates  of  state  3.  These  results  are  shown  in 
Figures  17  through  20.  Figures  17  and  18  present  the 
normalized  state  3 estimates  for  the  two  filters,  and 
Figures  19  and  20  show  the  normalized  percentage  error  in 
the  estimates.  The  extended  Kalman  filter  provided  a 
better  estimate  of  state  than  the  adaptive  filter  except 
during  the  initial  sample  period.  During  the  initial 
period,  the  adaptive  filter  produced  a more  accurate  esti- 
mate. This  result  is  more  significant  if  the  error  in  the 
initial  estimate  of  state  3 is  larger. 

In  Figures  21  through  33,  the  initial  error  in  the 
estimate  of  state  3 is  increased  by  setting  the  state  3 
value  equal  to  zero  for  the  two  filters.  During  the  initial 
period,  the  adaptive  filter  generated  more  accurate  esti- 
mates for  both  the  uncertain  parameter  and  the  two  system 
states.  This  is  due  to  the  fact  that  the  adaptive  filter 
operates  as  a fixed-time-lag  data  smoother  (i.e.,  30  data 
samples  are  acquired  prior  to  establishing  the  first  state 
estimate)  while  the  extended  Kalman  filter  operates  as  a 
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true  filter  (i.e.,  generating  the  first  state  estimate  as 
the  first  data  sample  is  obtained) . 

The  extended  Kalman  filter,  however,  provided 
slightly  better  estimates  during  later  periods  than  the 
adaptive  filter  and  did  not  display  any  signs  of  divergence 
(i.e.,  instability)  that  exist  for  the  extended  Kalman 
filter  in  some  problems.  For  example,  divergence  of  the 
extended  Kalman  filter  has  been  particularly  acute  in  orbit 
determination  problems  due,  often,  to  system  modeling  errors 
(see  Wolf  1968) . The  divergence  in  the  filter  can  occur 
when  the  error  covariance  matrix  becomes  very  small.  Since 
the  filter  gain  is  related  to  the  error  covariance,  it  also 
becomes  very  small  and  the  filter  estimate  becomes  de- 
coupled from  the  observational  sequence  such  that  the  esti- 
mate is  not  affected  by  a growing  observational  error. 

The  computational  time  for  operation  of  the  Kalman 
filter  on  the  Control  Data  Corporation  7600  computer  for 
the  300  measurement  samples  was  1.86  seconds.  The  times 
for  the  modified  Kalman,  extended  Kalman,  and  adaptive 
filters  were  1.90,  2.34,  and  9.56  seconds,  respectively. 

The  modified  Kalman  filter  provided  more  accurate  state 
estimates  in  about  the  same  time  as  the  Kalman  filter.  The 
extended  Kalman  filter,  however,  produced  more  accurate 


( 
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estimates  than  the  modified  Kalman  filter  with  a slight 
increase  in  computational  time.  The  adaptive  filter  did 
not  improve  the  accuracy  relative  to  the  extended  Kalman 
filter  with  the  exception  of  the  initial  sample  period  and 
consumed  a relatively  larger  computational  time.  The 
larger  time  for  the  adaptive  filter  was  due  to  the  number 
of  iterations  required  to  solve  the  nonlinear  estimation 
equation  for  state  3.  It  is  interesting  to  note,  however, 

j 

that  a significant  portion  of  computer  time  (i.e.,  about 
4 additional  seconds  for  each  filter)  was  required  just  for 
the  printing  of  results.  Therefore,  the  additional  compu- 
tational time  may  not  be  significant  compared  to  the  total 
computer  time. 

5.7  EFFECT  OF  NOISE  LEVEL 

The  effect  of  the  level  of  measurement  noise  on 
the  various  filter  estimates  is  illustrated  in  Figures  33 
through  72.  In  Figures  33  through  52  the  measurement 
noise  covariance  is  reduced  by  a factor  of  2,  and  in 
Figures  53  through  72  the  covariance  is  increased  by  a 
factor  of  2.  In  general,  the  error  in  all  the  filter 
estimates  was  reduced  when  the  noise  was  decreased  and  was 
increased  when  the  noise  was  increased.  However,  this 
reduction  and  enlargement  in  the  filter  error  was  more 
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significant  in  the  Kalman  and  modified  Kalman  filters.  The 
extended  Kalman  and  adaptive  filters  showed  little  dif- 
ference in  the  estimate  error  with  variation  in  measure- 
ment noise  level. 

5.8  DATA  ABERRATIONS 

The  effects  of  a gap  in  the  measurement  data  and 
an  erroneous  data  spike  were  examined  for  the  modified 
Kalman  filter  and  the  adaptive  filter.  Data  for  samples 
100  through  104  was  eliminated  from  the  filtering  process. 
The  effect  on  the  two  filter  estimates  is  shown  in  Figures 
73  through  82.  In  general,  the  error  in  the  estimate 
generated  by  the  modified  Kalman  filter  was  increased  by 
the  data  gap  beginning  at  the  time  of  the  first  missing 
sample.  The  increase  in  the  filter  error,  however,  was 
only  temporary  and  in  this  example  the  error  magnitude 
became  equivalent  to  previous  results  after  approximately 
20  samples.  In  contrast,  the  effect  of  the  data  gap  on 
the  adaptive  filter  was  small  and  differed  little  from 
results  generated  with  the  full  complement  of  data. 

The  data  at  sample  100  was  multiplied  by  a factor 
of  10  to  create  an  erroneous  spike  in  the  data  with  a 
variation  equivalent  to  more  than  100  times  the  standard 
deviation  of  the  measurement  noise.  Figures  83  through  92 


illustrate  the  effect  on  the  estimation  process  for  the 

two  filters.  The  error  at  sample  100  of  the  modified 

Kalman  filter  estimate  was  significantly  increased  and 

remained  greater  than  previous  results  for  approximately 

the  next  25  samples.  After  this  sample  period,  the 

% 

estimates  generated  by  the  modified  Kalman  filter  did  not 
differ  from  previous  results. 

Because  of  the  iterative  nature  of  the  adaptive 
filter,  the  jrror  for  sample  estimates  prior  to  sample  100 
was  increased  by  the  data  spike  at  sample  100.  In  addition, 
the  estimate  error  was  abnormally  higher  for  approximately 
the  next  50  to  60  samples.  After  this  period,  the  estimate 
error  was  equivalent  to  results  generated  without  the  data 
spike.  Due  to  the  dramatic  effect  of  the  data  spike  on  the 
estimate  error,  practical  use  of  the  adaptive  filter  may 
require  a pre-filter  to  limit  data  spikes  and  ameliorate 
the  effects  on  the  filter  estimation. 

The  data  spike  might  be  considered  as  noise  and, 
as  such,  be  indicative  of  a larger  noise  covariance.  A 
larger  covariance  could  then  be  substituted  into  the  filter 
algorithms.  The  effect  of  this  substitution  would  be  to 
decrease  the  gain  term  and  make  the  filter  less  responsive 
to  data  spikes,  and  thus,  provide  a better  estimate  for 
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these  samples.  The  filter,  however,  would  also  be  less 
responsive  to  samples  without  the  data  spike.  Since 
samples  without  data  spikes  reflect  a lower  noise  covariance 
and  because  the  filter  is  only  optimum  on  the  basis  of  an 
assumed  white  noise  process  with  known  covariance,  the 
filter  would  not  provide  the  best  estimates  for  these 
samples.  By  retaining  the  lower  noise  covariance  repre- 
sentive  of  the  measurement  noise  without  data  spikes,  better 
estimates  can  be  made  for  these  estimates  while  degraded 
estimates  will  be  made  of  samples  with  data  spikes.  Also, 


since  the  effect  of  the  data  spike  was  shown  in  this  ex- 
ample to  attentuate  with  time,  retaining  the  lower  noise 
covariance  may  be  more  desirable. 
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Initial  State  - 0.66  rad/sec  Sample  Rate  - 0.1  sec 


Figure  1.  Kalman  Filter  Estimate  of  State  1 Determined 
with  a Parameter  Error  in  the  System  Dynamics 
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Figure  4.  Adaptive  Filter  Estimate  of  State  1 Determined 
with  a Parameter  Error  in  the  System  Dynamics. 


+ 

t 

c 

1 c 

0 

•H 

I -w 

p 

p 

m 

i ifl 

•H 

•H 

> 

> 

V 

0) 

Q 

Q 

|: 

73 

Tj 

P 

R P 

(0 

IT) 

73 

I -a 

c 

s 

1 " 

(0 

P 

CO 

CO 

<u 

0) 

P 

f;  p 

io 

<0 

E 

CM 

£ 

•H 

. - 

•3 

p 

0 

4J 

03 

0) 

to 

H 

CM  CN 

03 

W 

, — 

\ 

(0 

u 

u 

'O 

0] 

3 

<d 

cu 

(0 

3 

C 

0) 

03 

p 

pH 

»H 

\ 

\ 

cu 

2 

73 

73 

(0 

<o 

in 

P 

p 

o 

a 

— ■ 

o 

a) 

o 

03 

rH 

pH 

• 

\ o 

o 

o 

o 

73 

o 

o 

0) 

(0 

o 

o 

1 

(0 

P 

• 

o 

u \ 

o 

• 

<D 

d) 

73 

v£> 

o 

0 

03 

3 

vD 

1 

c 

\ 

'O 

P 

• 

1 

(0 

O 

0 

•H 

n> 

m 

0 

a) 

p 

La 

lO 

1 

c 

0 

(0 

• 

<0 

c 

> 

0 

vO 

O 

P 

•H 

(0 

0 

a> 

vO 

c 

p 

•H 

u 

03 

• 

1 

3 

id 

P 

i 

o 

P 

> 

10 

43 

rH 

<u 

03 

0 

> 

03 

• 

1 

P 

c 

u 

o 

•H 

o 

<o 

0 

u 

0 

p 

0) 

E 

u 

a) 

z 

i 

0 

P 

•H 

p 

a> 

p 

n) 

P 

t3 

(0 

03 

p 

<u 

p 

P 

03 

c 

p 

•H 

c 

p 

w 

CO 

w 

•pH 

CO 

0 

a) 

Ifl 

N 

z 

E 

CO 

<v 

rH 

rH 

*H 

pH 

4) 

p 

<0 

3 

pH 

10 

E 

P 

0) 

Ifl 

•H 

•H 

m 

•H 

a) 

3 

rH 

E 

P 

■H 

c 

P 

E 

p 

P 

03 

a« 

•H 

•H 

p 

•H 

03 

(0 

P 

p 

c 

0 

c 

>i 

a> 

!3 

03 

H 

H 

z 

H 

CO 

2 

w 

W 

iiiiSBIl 


IpirilHBnBHEaBEeBS 

SiMainisKgiggei 

aiiSBSiiii: 


iie| 


SBSi 


Eiin 
(GSaS 

ummmm 

■Si! 


558 


nil 


PHMiiaiRi 


ElasiiaS' 
|»fi»S| 


BSiiiscissBBKaWMMHBHH 
n^iE£SSSa»SBK^Sii6ai:seBEBBSSCEl&iifi| 
fe‘He*3agSESBBSB'«I381EiSSBSSSgSB6lSBBii| 

fpiNliilBliiiiSIlilPiiillliiil 
l^iijiIiiilili*i*lii;-iiiiiiiliiisKp 

—»sHBafera»aca«KEBftBfaWk5a«^BBgrl 


SAMPLE 

Figure  5.  Kalman  Filter  Estimate  Error  in  State  1 Determined 
with  a Parameter  Error  in  the  System  Dynamics. 
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Figure  7.  Extended  Kalman  Filter  Estimate  Error  in  State  1 Determined 
with  a Parameter  Error  in  the  System  Dynamics. 
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Figure  11.  Extended  Kalman  Filter  Estimate  of  State  2 Determined  with 
Parameter  Error  in  the  System  Dynamics. 
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Figure  13.  Kalman  Filter  Estimate  Error  in  State  2 Determined  with 
Parameter  Error  in  the  System  Dynamics. 
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Figure  14.  Modified  Kalman  Filter  Estimate  Error  in  State  2 Determined  with 
Parameter  Error  in  the  System  Dynamics. 
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Figure  18.  Adaptive  Filter  Estimate  of  the  Uncertain 
Parameter  in  the  System  Dynamics. 
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Figure  21.  Extended  Kalman  Filter  Estimate  of  State  1 Determined  from  a Zero  Initial 
Estimate  for  the  Uncertain  Parameter  in  the  System  Dynamics. 
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Figure  27.  Extended  Kalman  Filter  Estimate  Error  in  State  2 Determined  from  a Zero 
Initial  Estimate  for  the  Uncertain  Parameter  in  the  System  Dynamics. 
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Figure  28.  Adaptive  Filter  Estimate  Error  in  State  2 Determined  from  a Zero 

Initial  Estimate  for  the  Uncertain  Parameter  in  the  System  Dynamics 
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Figure  33.  Kalman  Filter  Estimate  of  State  1 Determined 
with  Reduced  Measurement  Noise. 
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Figure  34.  Modified  Kalman  Filter  Estimate  of  State  1 Determined 
with  Reduced  Measurement  Noise. 
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Figure  43.  Extended  Kalman  Filter  of  State  2 Determined 
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Figure  44.  Adaptive  Filter  Estimate  of  State  2 Determined 
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Modified  Kalman  Filter  Estimate  Error  in  State 
Determined  with  Reduced  Measurement  Noise. 
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Determined  with  Reduced  Measurement  Noise. 


Parameter  Value  - 7.4  rad/sec  Estimate  Error  

Initial  Value  - 6.3  rad/sec  Plus  Estimate  Standard  Deviation  • 

Normalizing  Constant  - 7.4  rad/sec  Minus  Estimate  Standard  Deviation 

Measurement  Noise  Covariance  - 0.00025  (rad/sec)2 
Sample  Rate  - 0.1  sec 


Figure  51.  Extended  Kalman  Filter  Estimate  Error  in  the  Uncertain  Parameter  in  the 
System  Dynamics  Determined  with  Reduced  Measurement  Noise. 


Figure  52.  Adaptive  Filter  Estimate  Error  in  the  Uncertain  Parameter  in  the  System 
Dynamics  Determined  with  Reduced  Measurement  Noise. 
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Initial  Estimate  - 0.63  rad/sec  True  State  


Fil 
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Figure  55.  Extended  Kalman  Filter  Estimate  of  State 
Determined  with  Increased  Measurement  No 


Figure  56.  Adaptive  Filter  Estimate  of  State  1 Determined 
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Figure  57.  Kalman  Filter  Estimate  Error  in  State  1 Determined 


Initial  State  - 0.66  rad/sec  Sample 
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Figure  58.  Modified  Kalman  Filter  Estimate  Error  in  State 
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Figure  63.  Extended  Kalman  Filter  Estimate  of  State  2 Determined 
with  Increased  Measurement  Noise. 
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Figure  69.  Extended  Kalman  Filter  Estimate  of  the  Uncertain  Parameter  in  the  System 
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Figure  71.  Extended  Kalman  Filter  Estimate  Error  in  the  Uncertain  Parameter  in  the 
System  Dynamics  Determined  with  Increased  Measurement  Noise. 
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Figure  72.  Adaptive  Filter  Estimate  Error  in  the  Uncertain  Parameter  in  the  System 
Dynamics  Determined  with  Increased  Measurement  Noise. 
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Figure  81.  Adaptive  Filter  Estimate  of  the  Uncertain  Parameter  in  the  System 
Dynamics  Determined  with  Missing  Measurement  Data. 
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Figure  82.  Adaptive  Filter  Estimate  Error  in  the  Uncertain  Parameter  in  the  System 
Dynamics  Determined  with  Missing  Measurement  Data. 


Initial  State  - 0.66  rad/sec  Measurement  Data  at  Sample  100  Multiplied  by  10 
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Figure  84.  Adaptive  Filter  Estimate  of  State  1 Determined 
with  an  Erroneous  Measurement  Data  Spike. 
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Figure  88.  Adaptive  Filter  Estimate  of  State  2 Determined 
with  an  Erroneous  Measurement  Data  Spike. 


Figure  89.  Modified  Kalman  Filter  Estimate  Error  in  State  2 

Determined  with  an  Erroneous  Measurement  Data  Spike 


CHAPTER  6 


L 


CONCLUSION 

This  chapter  provides  a summary  of  the  research 
results  and  recommendations  for  future  investigations. 

6.1  SUMMARY  OF  RESULTS 

In  the  dissertation,  two  filter  techniques  are 
proposed  for  estimation  of  the  states  of  a linear  dynamic 
system  using  noise-corrupted  measurements  when  parameters 
in  the  system  model  are  not  exactly  known.  The  first 
technique  is  a modification  of  the  Kalman  filter.  The 
method  incorporates  the  parameter  uncertainties  in  the 
filter  by  adding  additional  terms  to  the  Kalman  gain  equa- 
tion. In  general,  when  the  parameter  uncertainties  exist 
in  the  state  equation,  the  effect  of  these  terms  is  to 
increase  the  filter  gain  since  the  measurement  is  a more 
accurate  assessment  of  the  state  than  the  predicted  state. 
Conversely,  when  the  parameter  uncertainties  exist  in  the 
measurement  equation,  the  gain  is  decreased  since  the 
predicted  state  is  a more  accurate  appraisal  relative  to 
the  measurement. 

The  second  technique  is  an  adaptive  method  based  on 
the  properties  of  the  innovation  sequence  of  the  optimum 
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Kalman  filter.  This  method  makes  a maximum  likelihood 


Bayesian  estimate  of  the  uncertain  parameters  operating 
iteratively  in  conjunction  with  a Kalman  filter.  The 
Kalman  filter  uses  the  current  estimate  of  the  parameters 
to  generate  the  innovations  and  other  information  utilized 
by  the  Bayesian  estimator  of  the  parameters.  With  the 
improved  values  for  the  parameters,  the  Kalman -filter  is 
again  operated  to  regenerate  the  state  estimates. 

Numerical  results  for  an  example  problem  were 
generated  for  the  modified  Kalman  filter  and  the  adaptive 
filter.  For  comparison,  data  was  also  generated  for  the 
Kalman  and  extended  Kalman  filters.  In  general,  the 
modified  Kalman  filter  produced  a more  accurate  estimate 
than  the  Kalman  filter  in  about  the  same  computational  time. 
The  extended  Kalman,  however,  provided  a significantly 
better  estimate  than  either  the  Kalman  or  modified  Kalman 
filter  by  estimating  the  uncertain  parameter  but  did  con- 
sume slightly  more  computational  time.  The  accuracy  of 
the  adaptive  filter  was  slightly  less  than  that  of  the 
extended  Kalman  filter  except  during  the  initial  sample 
period  and  the  computational  time  was  significantly  greater. 
Since  the  computational  time  is  only  a few  seconds,  the 
larger  computational  time  may  not  be  a disadvantage  for 
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the  adaptive  filter  in  some  problems.  In  addition,  the 
adaptive  filter  may  be  more  useful  since  the  extended 


Kalman  filter  is  sometimes  plagued  by  divergence.  The 
basic  reason  that  the  estimate  diverges  from  the  actual 
state  is  that  the  gain  in  the  filter  algorithm  approaches 
zero  too  rapidly.  The  multiplication  of  the  gain  and  the 
innovation  (i.e.,  the  measurement  residual)  in  the  filter 
algorithm  thus  has  a negl igible  effect  on  the  estimate  even 
though  the  innovation  may  be  quite  large.  Hence,  the 
estimate  becomes  decoupled  from  the  observation  sequence 
and  is  not  affected  by  the  increasing  measurement  residuals. 
On  the  other  hand,  by  designing  the  adaptive  filter  to 
provide  a maximum  likelihood  estimate,  the  measurement 
residuals  form  the  basis  for  the  estimate.  Operation  of 
the  filter  attempts  to  minimize  the  measurement  residuals. 
Consequently,  the  filter  attempts  to  adjust  and  correct  the 
estimates  with  increasing  measurement  residuals  and  thus 
ensure  filter  stability.  The  proof  of  this  behavior  is  not 
considered  here,  but  could  be  addressed  in  future  studies. 

When  the  level  of  measurement  noise  was  increased 
or  decreased,  the  errors  in  the  estimates  generated  by  the 
four  filters  tended  to  increase  or  decrease  correspondingly. 
However,  the  effect  was  more  significant  on  the  results  of 
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the  Kalman  and  modified  Kalman  filters.  The  extended 
Kalman  filter  and  adaptive  filter  showed  relatively  little 
difference  with  the  variation  in  measurement  noise  level. 

The  effect  of  missing  data  on  the  modified  Kalman 
filter  was  analyzed  and  found  to  increase  the  estimation 
error.  The  increase,  however,  was  only  temporary  in  the 
example  and  became  equivalent  to  results  generated  with  the 
complete  set  of  data  after  approximately  20  samples.  The 
effect  of  the  data  gap  on  the  adaptive  filter  was  relatively 
small  and  differed  little  from  results  generated  with  the 
full  complement  of  data. 

An  erroneous  data  spike  was  also  analyzed  and 
produced  a significant  increase  in  the  estimation  error  of 
the  modified  Kalman  filter,  which  subsequently  affected 
approximately  25  state  estimates.  After  this  period,  the 
estimates  returned  to  the  values  generated  without  the 
data  spike.  Because  of  the  iterative  nature  of  the 
adaptive  filter,  the  error  in  the  adaptive  filter  esti- 
mates prior  to  the  measurement  data  spike  was  increased 
and  the  estimation  error  was  abnormally  higher  for  50  to 
60  estimates  after  the  data  spike.  The  estimates  did, 
however,  return  to  the  values  generated  without  the  data 
spike.  Due  to  the  dramatic  effect  of  the  data  spike  on  the 


225 


estimation  error,  practical  use  of  the  adaptive  filter  may 
require  a pre-filter  to  limit  the  effects  of  data  spikes. 
Additional  experimental  or  data  analyses,  however,  would 
be  required  to  develop  and  establish  the  design  of  the 
pre-filter  or  techniques  which  treat  effectively  data 
spikes . 

Tables  1 through  4 summarize  the  numerical  results. 
In  general,  the  Kalman  filter  can  be  applied  to  the  problem 
of  state  estimation  for  systems  with  uncertain  parameters 
provided  that  the  estimates  meet  the  desired  accuracy 
requirements.  If  greater  accuracy  is  required,  then  the 
modified  Kalman,  the  extended  Kalman  or  the  adaptive  filter 
might  be  considered.  The  desirability  of  each  filter, 
however,  must  be  addressed  on  a case  by  case  basis.  The 
extended  Kalman  filter  can  produce  greater  accuracy  than 
the  modified  Kalman  filter  provided  that  divergence  of  the 
extended  Kalman  filter  estimates  is  avoided.  The  adaptive 
filter  can  produce  about  the  same  accuracy;  however,  the 
computational  time  of  the  filter  may  be  greater. 

6.2  RECOMMENDATIONS  FOR  FUTURE  INVESTIGATION 

There  are  several  areas  of  investigation  that 
could  be  considered  for  future  effort.  These  topics  are 
described  in  the  following  paragraphs. 


r 
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First,  the  innovation  sequence  of  the  Kalman  filter 
was  used  in  the  adaptive  filter  to  generate  state  esti- 
mations based  on  a predefined  schedule  for  the  measure- 
ments. If  the  frequency  of  measurements  is  increased  or 
decreased,  the  magnitude  of  values  of  the  innovation 
sequence  and  its  accompanying  covariance  would  be  expected 
to  increase  or  decrease  correspondingly.  This  leads  to  a 
change  in  the  accuracy  of  the  state  estimates.  Since  the 
taking  and  processing  of  measurements  represent  a certain 
effort  or  cost,  a tradeoff  exists  between  the  estimate 
accuracy  and  the  cost  of  taking  the  measurements.  Thus, 
by  utilizing  the  innovation  sequence,  an  adaptive  method 
might  be  developed  that  determines  the  measurement  schedule 
while  maintaining  a certain  degree  of  estimation  accuracy. 

Secondly,  the  adaptive  filter  consumed  a relatively 
larger  amount  of  computational  time  than  other  filters 
analyzed.  Therefore,  methods  for  reducing  the  computa- 
tional time  could  be  investigated.  For  example,  certain 
quantities  might  be  precomputed,  or  summetry  in  matrix 
relationships  might  be  exploited. 

. Third,  the  divergence  problem,  or  instability,  of 
the  extended  Kalman  filter  was  noted  in  Chapter  5. 

Although  some  studies  of  stability  have  been  performed 


for  the  Kalman  filter,  studies  of  the  stability  of  the 
extended  Kalman  filter  are  almost  nonexistent.  Therefore, 


study  efforts  on  the  extended  Kalman  filter  stability 
should  be  performed  to  guide  the  determination  of  appro- 
priate applications  for  the  extended  Kalman  filter.  In 
addition,  stability  of  the  adaptive  filter  utilizing  the 
maximum  likelihood  estimator  could  also  be  addressed  in 
such  a study. 

Finally,  the  simultaneous  problem  of  optimum  con- 
trol and  state  estimation  for  systems  with  certain  para- 
meters might  be  considered.  In  this  study,  only  the 
problem  of  estimation  was  addressed.  Probably  the  easiest 
as  well  as  an  important  case  of  optimal  control  to  be 
examined  would  be  linear  feedback  control.  Another  area 
of  potential  study  is  the  development  of  a separation 
theorem  between  the  optimal  control  and  the  estimate  for 
systems  with  uncertain  parameters.  This  study  might 
initially  consider  the  report  by  Snyder  (1977) . 
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APPENDIX  A 

DISCRETE  REPRESENTATION 
OF  CONTINUOUS  DYNAMIC  SYSTEMS 

This  appendix  will  develop  the  discrete  system 
model  that  is  equivalent  to  the  continuous -time  descrip- 
tion. Assume  that  the  continuous -time  model  of  the  system 
is  given  by  the  stochastic  differential  equation 

dx  = Fx  dt  + Bu  dt  + G d/3  (A-l) 

or  the  equivalent,  through  mathematically  less  precise, 
linear  differential  equation 

x = Fx  + Bu  + Gw  (A-2) 


where  w is  white  Gaussian  noise  and  formally  the  derivative 
of  Brownian  motion  [3 . The  white  noise  is  a zero-mean 
process  with  covariance  Q(t)  8(t-r)  where  Q(t)  is  chosen 
to  duplicate  the  low  frequency  power  spectral  density  of 
the  actual  noise  entering  the  system.  The  statistics  of 


the  noise  are  expressed  by  the  equations 


E[w(t1)wT(t2)]  = Q(tL)  8(t1-t2) 


(A-4) 


It  is  assumed  that  the  measurements  are  taken  at 
discrete  instants  and  are  of  the  form 


z . = H.X.+  v. 

i ill 


(A- 5) 


Furthermore,  it  is  assumed  that  a digital  computer  will 
provide  the  control  input,  so  that  u(t)  will  be  piecewise 
constant.  That  is,  a measurement  would  be  taken,  the 
information  processed,  and  a control  input  created  and 
held  constant  until  the  following  sample  time. 

The  solution  to  Equation  (A-l)  or  (A-2)  is 

x(t)  = ^(t,t^)x(t^)  + J 0(t,r)B(r)u(T)dt 


+ f 0(t,  r)G(r)dw(r) 


(A-6) 


where  the  state  transition  matrix  0(t,t^)  satisfies  the 
differential  equation  and  initial  condition 


0(t,ti)  = F(t)0(t,ti) 


(A-7) 
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I 


(A-8) 


Since  the  control  is  held  constant  over  a given  sample 
period,  the  solution  within  a single  sample  period  is 

X ( t)  = 0(t,ti)x(ti)  + i|»(t,ti)u(ti) 


/*<t. 


r)  G (r)  dw  (r) 


(A-9) 


rti 


where  the  matrix  <^(t,t^)  is  defined  as 


/*  ^ 

* J 0(t,r)B(r)dr 


(A-10) 


Differentiating  this  with  respect  to  time  yields  the 
equivalent  defining  relations 


</,(t,ti)  = B (t)  + F(t)^(t,ti) 


/-ti 

^(ti,ti)  = J 0(t,T)B(T)dr 


Thus,  if  an  equation  of  the  form 


= 0 


(A-ll) 


(A-12) 


xi+i  ■ *i+i/ixi + Wiui  * ri+i/iwi  (R-13) 


is  to  duplicate  the  state  at  the  (i+1)  sample  time,  x(ti+1) , 
as  given  by  Equation  (A-9) , then  it  can  be  seen  that 


0i+l/i  ^tti+l»ti)  (A-14) 

‘/'i+l/i  = ^(ti+i,ti)  (A-15) 


For  filter  applications,  the  last  term  in  (A-13) 
need  not  be  evaluated.  Rather,  an  expression  for  the 
covariance  of  its  contribution  is  required: 


fci+l 


WA'Wi  -ft.  *(ti+i'r>G'r>°'r>GT(T>*T<ti+rT>dT 

(A- 16) 


or 


Ji+l/iQi/i+l/i'’  W^ti+l,ti^  (A-17) 

where  Equation  (A-17)  serves  as  a definition  of  W(ti+^,t^). 
As  in  the  previous  case,  a more  convenient  form  for 
evaluation  is  obtained  by  differentiating  Equation  (A-16) 
to  yield 
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Thus,  to  obtain  the  equivalent  discrete  model  for 


a continuous -time  system.  Equations  (A-7) , (A-ll) , and 

(A-18)  are  integrated  forward  from  the  initial  conditions 

(A-8) , (A-12) , and  (A-19)  to  time  ti+1.  The  required 

T 

matrices  </'i+i/i*  and  /i+i/iOi/i+i/i  are  then 

determined  from  Equations  (A-14) , (A-15) , and  (A-17) , 

respectively. 

The  influence  of  uncertain  parameters  in  F(t)  or 

B (t)  upon  the  discrete  model  can  be  expressed  analytically 

in  simple  cases.  In  more  complex  situations,  the  dependence 

can  be  found  by  numerical  integration  of  these  equations 

for  various  values  of  the  parameters,  from  which  functional 

relationships  between  the  elements  of  *1+1/1 • and 

T 

^i+l/iQ/i+l/i  201,3  t3ie  Paran»eter  values  can  be  established. 
If  the  system  is  assumed  to  be  time  invariant  and  noise 
inputs  are  assumed  to  be  stationary  over  N samples,  a 
single  set  of  integrations  suffices  for  all  sample  periods. 


237 


238 


1 


APPENDIX  B 


MINIMUM  ERROR  VARIANCE  ESTIMATE 


Consider  a filter  estimate  x^  that  is  to  be  defined 
from  measurements  z^,  z^.-.z^  so  that  it  minimizes  the 
trace  of  the  error  variance  of  the  estimate.  Equationally, 
this  is 


tr  ^[(xj^-x*)  (xi-x*)T]j  = minimum  (B-l) 


or 


E [(x^-x*) T (xi-x*)]  = minimum 


(B-2) 


To  derive  this  estimate,  the  error  variance  must  first  be 
written  in  terms  of  the  density  function 

E[(xi-x*)T(xi-x*)]  = //  ./  (Xi-x^K-x*) 


p(x.,z....z  ) dx.dz,...dz. 
* i 1 i i 1 i 


(B— 3 ) 


where  each  integral  sign  represents  a multiple  integral 
since  x^,  z^...z^  may  be  vectors.  The  density  function  cam 
be  written  using  Bayes'  theorem  as 

239 


E I x . I z 


By  definition,  this  quantity  is  positive.  Therefore,  to 
minimize  the  multiple  integral,  it  is  sufficient  to  choose 
x*  to  minimize  the  integrand  given  by  Equation  (B-S) . Only 


the  last  term  involves  x^  and  the  smallest  value  it  can 
assume  is  zero.  Thus,  the  minimizing  estimate  is  given  by 


x*  = EOJZ-L...ZJ 


(B-7) 


The  Kalman  filter  estimate  x^  for  linear  Gaussian  systems 
is  also  given  by 


*i  = EC*if  zx*  • 


(B-8) 


Therefore,  the  Kalman  filter  estimate  is  a minimum  error 
variance  estimator  (i.e.,  = x*) . 
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APPENDIX  C 


GAIN  RELATIONSHIP 
FOR  PERFECT  MEASUREMENTS 

If  the  measurements  are  directly  related  to  the 
state  in  such  a manner  that  aHi  is  zero,  the  modified 
filter  gain  is  given  by  the  expression 

Gi  = ^(11^  + Ri)"1  (C-l) 

When  the  noise  covariance  R^  is  zero,  the  measurements 
determine  the  state  exactly  and  the  filter  gain  becomes 

Gi  - (C-2) 

or 


HiGi  = 


-1 


(C-3) 


= I 


Premultiplying  by  hT  and  performing  the  inverse  operation 
gives 
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APPENDIX  D 


DERIVED  DENSITY  THEOREM 


Let  z and  V be  random  N vectors  where  z equals  f(v) 
and  v equals  f-1(z)  with  f and  f-1  continuously  differen- 
tiable. For  an  arbitrary  set  S 


Pr(zfS)  = Pr 


(f  iv)es) 


= Pr  (t'€f*1(S)) 


(D-l) 


where  Pr  ( ) is  the  probability  of  ( ) . In  terms  of  the 
probability  density  function  of  z,  pz (z) , this  probability 
can  be  written 


Pr(z€S)  = / pz(z)dz 


(D-2) 


In  addition,  the  right  half  of  Equation  (D-l)  can  be 
written  using  the  probability  density  function  £v(v)  as 


Pr[»'€f~1(S)]  = f p„  (v)dv 


f _1(s) 

,/P„|£-1(z,|||^l||dz 


(D-3) 


f [pzu>  - p„(f"1U>)ll  2CilSl||]  dy  - 0 (D-4) 


Since  this  equation  must  hold  for  any  arbitrary  region  of 
integration  no  matter  how  small  or  large,  the  integrand 
must  be  zero.  Therefore,  the  densiLy  function  of  z 
derived  from  the  density  function  of  v is 


Pz(z)  = *v  f 


■1(z)|ll 


df-1 (z) 
dz 


(D-5) 
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APPENDIX  E 


r 


EVALUATION  OF  FISHER  INFORMATION  MATRIX 

The  Fisher  information  matrix  developed  in 
Chapter  4 is  given  by  the  expression 


N 


F(<V“/> 


4E 


£1  v-l  ±± 


+ E 


T -1  5Vi  -1  T -1  aVi  -1 

-L  r\  1 1 x 1 ~ 11 

a°k  dai 


trfv'1  £!i)  trfv'1  £l> 

8“k 


(E-l) 


Determination  of  F (a^,a^)  requires  an  evaluation  of  the 
second  term  of  this  expression.  This  evaluation  is  accom- 
plished by  letting  the  innovation  be  defined  in  terms 
of  a vector  A whose  components  are  independent  Gaussian 
variables  with  zero  mean  and  unity  variance.  The  desired 
relationship  is 


"i  - A 


(E-2) 


so  that 
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r 


E^^i)  = /w~  E[AAT]  >/vf 


= /v“  i/v~ 


= V, 


(E-3) 


which  is  consistent  with  the  defined  variance  of  v . Also, 

i 

let  symmetric  matrices  a S'  and  a,S"  be  defined  by 

3c  jL 


„ S'  = JT?  VTlf^VTl  JT. 
ak  v i i i v i 


— T _2.  / 

S"  = yvT  v±L  /v^ 

0t£  1 Oa £ 1 1 


(E-4) 


(E-5) 


The  second  term  of  the  information  matrix  becomes 


t -i  dvL  _i  T _i  av±  _i 
V.v.  — - v.  v.v.v.  — - v.  v. 
1 1 dc^  1 1 1 1 daL  1 * 


= E 


= E 


AT  S'AAT  s"a 
. °k  a/ 


tr  ( S*AAT)tr(  S"AAT) 
“k  a£ 


(E-6) 
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Since  the  components  of  A,  A^f  are  independent,  first  and 
third  moments  are  zero  and  second  and  fourth  moments  are 
1 and  3,  respectively.  Therefore,  the  expected  value 
above  can  be  written 


E tr  ( S 'AAT)  tr  ( S'AAT)  = tr  ( S * ) tr  ( S")  + 2tr(  S'  S") 

“k  “/  J ak  al  “k 

(E-7) 


The  second  term  of  the  information  matrix  then  becomes 


T -1  dVi  _i  t -1  dVi  -1 

E v.V.  i.  V.  V.  V.V.  V.  V.  = 

1 1 da^  1 1 1 1 daL  1 1 


tr  (/v*  V-1  VT1  TvT)  tr(  v/v*  vT1  fli  v’1  V^T ) 


2tr  (VvT  VTl  fZi  VT1  Jv?  V,"1  fXi  VT1  ) 

'll  1 111  1 1 / 


= trfvT1  fliWV1  !li)  + 2tr(vT1  fZi  vT1  ) (E-8) 

' 1 A.  / ' 1 / ' 1 1 / 


da  da/ 

k x 


Substituting  into  Equation  (E-l)  and  utilizing  the  sample 
data  yields 
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which  completes  the  evaluation  of  the  Fisher  information 
matrix. 
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APPENDIX  F 


PARAMETRIC  LIMITS 
OF  SPINNING  BODIES 

The  basic  parameters  of  spinning  bodies  are:  the 

angular  momentum  of  the  body  M;  the  kinetic  energy  of  the 

body  T;  and  the  moments  of  inertia  for  the  three  principal 

axes  of  the  satellite  1^,  I2»  and  I^.  The  moments  of 

inertia  define  body  shape  constraints  while  the  angular 

momentum,  kinetic  energy  and  moments  of  inertia  dictate 

limits  of  rigid  body  motion. 

F.l  MOMENT  OF  INERTIA  CONSTRAINTS 

The  moments  of  inertia  about  each  axis  are  the 

integrated  products  over  the  body  of  the  mass  of  each 

element  in  the  body  and  the  square  of  the  distance  from 

the  axis  to  the  mass  element.  The  minimum  moment  of 

inertia  1^  is  about  the  b^  axis,  the  intermediate  value 

I2  is  about  the  b2  axis,  and  the  maximum  value  I3  is 

about  the  £3  axis  (see  Figure  F-l) . If  r^  is  the  generic 

a , 

scalar  component  along  the  b^  axis  to  a mass  element  dm, 

A 

and  r2  and  r^  are  generic  scalar  components  along  the  b2 
and  63  axes  so  that  the  generic  vector  is  r_  = r^b^  + 
r2b2  + *2—3'  t^ien  moments  of  inertia  are 
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rX  =/<r 2 + rl>  ^ 


(F-l) 


z2  mfizi  * rl>  *" 

I3  =y (r^  + rj)  dm 


Adding  Equations  (F-l)  and  (F-2)  gives 


1 + *2  =/(r  1 + r2>  <*"  + 2/r3  dm 


or 


h + *2  - X3  + 2/r 3 dm 


Equation  (F— 5 ) leads  to  the  inequality 


or 


I1  + I2  ” *3 


*3  *3 
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(F-2) 


(F-3) 


(F-4) 


(F-5) 


(F-6) 


(F— 7 ) 


■\r 


This  condition  is  illustrated  by  the  lined  area  in 


Figure  F-2.  The  equality  of  Equation  (F— 7 ) represents  the 
limiting  body  shape,  which  is  a flat  plate  with  negligible 
thickness . 


1.0 


Figure  F-2. 


0 


1.0 


Constraint  on  Moments  of  Inertia  (I^+^^Ig)  • 


The  moments  of  inertia  are  defined  as 


X3  * X2  * X1 


or 


(F-8) 


(F-9) 


— . ... 


This  condition  is  illustrated  by  the  lined  area  in 
Figure  F-3. 


0 13/13  1.0 


Figure  F-3.  Moment  of  Inertia  Constraint. 

Superimposing  Figure  F-3  on  Figure  F-2  gives 
Figure  F-4  which  shows  the  two  constraints  on  the  ratios 
13/13  and  I3/I3 • The  double-lined  area  in  Figure  F-4 
reveals  the  boundary  limits. 


The  ratio  1 3/1.3  b°unde(i  *>y  the  values  0.5  and  1.0,  while 
the  bound  for  I3/I3  depend  on  the  value  of  13/13.  The 
upper  extreme  is  13/13  or  a symmetrically  shaped  body 
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k ; 


and  the  lower  limit  is  1 - I2//l3  or  a plate  with 

negligible  thickness. 

F.2  RIGID  BODY  MOTION  CONSTRAINTS 

The  general  motion  of  a spinning  rigid  body  of 
symmetry  in  the  absence  of  applied  moments  can  be  described 
by  steady  precession  or  coning  of  the  symmetry  axis  (or 
spin  axis)  about  the  angular  momentum  vector  M,  which  is 
fixed  in  inertial  space  (see  Figure  F-5) . The  term  spin 
axis  is  applied  to  the  body  axis  with  which  the  angular 
velocity  vector  is  nominally  aligned  for  steady  spinning 
motion,  with  no  precession.  For  an  axisymmetric  body, 

the  symmetry  axis  must  be  selected  as  the  spin  axis . - 

Angular  Momentum  M 


Figure  F-5.  Vector  Representation  for  Free  Body  Motion. 

The  angle  0 between  the  spin  axis  and  the  angular 
momentum  vector  is  called  the  nutation  angle  or  coning 
angle  and  remains  constant  for  a rigid  body  of  symmetry 
during  the  precession  cycle.  An  asymmetric  body  also 
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exhibits  precession  of  the  spin  axis  but,  in  addition, 
displays  an  oscillation  in  $ during  the  precession  cycle 
that  is  called  nutation.  The  frequency  of  oscillation 
in  nutation  is  twice  that  of  precession. 

(The  above  definitions  for  precession  and  nutation 
are  based  on  etymology  and  classical  usage.  In  modern 
applications  of  space  sciences  and  gyrodynamics  other 
usages  are  common.  Unfortunately,  with  regard  to  axisym- 
metric  bodies  the  word  nutation  is  often  applied  to  the 
coning  motion  characterized  here -as  precession. ) 

A A A 

With  the  origin  of  body  axes  b^,  bg/  and  b^ 
coinciding  with  the  body  center  of  mass,  the  rotational 
equations  of  motion  for  a rigid  body  in  the  absence  of 
applied  torques  are 

1ld)l  ” (I2"I3>  0>20i2  = 0 (F-10) 

12<b2  “ (W  "lw3  = ° (F_11) 

13<i3  “ (Il"I2)  °il<U2  = ° (F"12) 


where  a>^,  (o^ 


, and  (d ^ are  the  angular  velocities  about  the 
respective  body  axes  and  I I and  are  the  principal 
moments  of  inertia  defined  so  that  I >1 --I.  • Multiplying 

M X 

Equations  (F-10)  through  (F-12)  by  o>1#  g>2,  and  o>3, 
respectively,  and  adding  and  then  integrating  yields 

11U>1  + 12OJ2  + 13UJ3  = 2T  (F-13) 

where  the  constant  of  integration  2T  is  twice  the  kinetic 
energy  of  rotation.  In  an  analogous  way,  multiplying 
Equations  (F-10)  through  (F-12)  by  12^2'  anc^  I3&,3' 

adding  and  integrating  gives 

(I^)2  + U2w2)2  + U3W3) 2 = M2  (F-14) 

where  M is  the  angular  momentum. 

With  the  equations,  the  motion  of  a rigid  body 
spinning  about  the  axis  of  maximum  inertia  can  be  stated 
by  the  conditional  equation 

2T  I3  > M2  £ 2T  I2  (F-15) 

or 
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(F— 16) 


1 > 


M" 


2T  I. 


The  geometric  interpretation  of  Equation  (F-16)  can  be 

illustrated  by  considering  the  nutation  angle.  The  upper 
2 

limit  of  M /2T  I3  implies  that  the  nutation  angle  is  zero 
(i.e.,  pure  spin  about  the  maximum  inertia  axis)  while  the 
lower  limit  implies  that  the  nutation  angle  approaches  a 
maximum  value  of  90  deg.  Figure  F-6  illustrates  the 

variation  in  the  limit  of  the  nutation  angle  for  values  of 

, 2 2 
I2/l2  and  M /2T  I^.  As  the  value  of  M /2T  decreases 

from  one,  the  nutation  angle  increases  from  0 to  90  deg. 

For  a symmetric  body  (i.e.,  1^  = I^)  with  constant 

spin  about  the  axis  of  maximum  moment  of  inertia  o»^. 

Equations  (F-10)  to  (F-12)  can  be  used  to  show  that  the 

spin  axis  moves  with  respect  to  the  body  axes  about  the 

angular  momentum  vector  at  a constant  nutation  angle  with 

a frequency  Si  given  by 


Q = 


(I3-I1)  w3 


(F-17) 
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In  addition,  since  the  moments  of  inertia  are  defined  so 


that  I3  is  greater  than  1^,  the  spin  axis  will  rotate 
about  the  fixed  angular  momentum  vector  in  retrograde 
precession. 
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